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CHAPTER I 


INTRODUCTION 

In the analysis of aircraft motion, wind effects are a key 
parameter, and thus deserve careful consideration. Past studies of 
this type have taken only constant winds into account. However, the 
analysis conducted in association with this contract utilizes a 
numerical tabulation of eleven thunderstorm outflows (converted to 
horizontal spatial coordinates with Taylor's hypothesis) measured with 
the 1,500-ft (500 m) instrumented, meteorological tower at the National 
Severe Storms Laboratory in Norman, Oklahoma [1], These tabulated 
data, when combined with two-dimensional equations of aircraft motion, 
allow for the computation of the aircraft behavior at any position in 
(x,z) plane within the flow regime. 

The necessity for a study of this nature precipitates from the 
alarming rate of increase in fatal and near-fatal accidents which have 
been attributed to wind shear. Presented below are several prime 
examples of wind shear related accidents: 

1. Eastern Airlines B-727 crashed on approach to 
John F. Kennedy International Airport Runway 22L on June 24, 1975, 
resulting in 113 fatalities. Probable cause: adverse winds associated 

with a very strong thunderstorm along the ILS (Instrument Landing 
System) glide path [2]. 


Numbers in brackets correspond to similarly numbered 
references in the Bibl iography . 



2. Continental Airlines B-727 crashed after takeoff on 

August 7, 1975, on Runway 35L at Denver-Stapleton International Air- 
port, seriously injuring 15 persons. Probable cause: aircraft 

encounter with severe wind shear (generated by outflow from a thunder- 
storm over the departure path) at an altitude and airspeed which 
precluded recovery to level flight [3]. 

3. Iberian Airlines DC-10-30 crashed on December 19, 1973, 

500 ft (152 m) short of Runway 33L at Boston-Logan International 
Airport, injuring 13 passengers. Probable cause: increased rate of 

descent induced by an encounter with a low-altitude wind shear and the 
captain's inability to compensate for it [4]. 

4. Delta Airlines DC-9 impacted a flood-control dike 785 ft 

(240 m) from the Runway 20 threshold at Chattanooga Municipal Airport 
on November 27, 1973, injuring 42 persons. Probable cause: excessive 

rate of descent initiated by a wind shear condition which existed in 
the lower levels of the approach path [5]. 

In numerous NTSB (National Transportation Safety Board) air- 
craft accident reports related to weather, the low-level wind shear 
hazard theme, particularly during takeoff and landing maneuvers below 
1,500 ft (500 m), recurs. In all likelihood, many more aircraft 
accidents attributable to wind shear have not been recorded as such, 
due to the relatively new awareness of this phenomenon. 

Realistic mathematical models of wind shear conditions are 

needed for use in flight simulators for the training of flight crews in 

aircraft handling techniques when encountering this type of situation. 

Fast time computer analysis to relate the potential hazards of wind 

shear to different types of aircraft and control systems is another 
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desired achievement for safety-design considerations. To do this, a 
computer code for investigating flight through thunderstorms has been 
developed which solves two-dimensional, non-linear equations of air- 
craft motion, incorporating with it the influence of wind shear. The 
equations of motion with wind shear are developed in Appendix A. Full 
documentation of the computer code is given in Appendix B. A further 
refinement of the main code adds randomly-generated turbulence effects 
to the existing variable wind profiles. Using this approach, a 
realistic simulation of aircraft motion while passing through thunder- 
storm outflows is determined, as the study will show. 

In obtaining the results presented, a careful study of air- 
craft flight paths through eleven modeled thunderstorm outflows was 
conducted. As a consequence of this analysis it was found that the 
effect of shear can be considered in terms of the aircraft's deviation 
from the glide path and touchdown point, and from the trajectory it 
follows. The trajectory, particularly with fixed controls, often 
traced the form of a phugoidal oscillation of varying amplitude, 
period, and frequency. The phugoidal oscillations are discussed 
further in Chapter IV. 

In all, 125 flights were simulated using the aerodynamic 
coefficients and physical data characteristic of the following air- 
craft: B-747, DC-8, augmentor-wing STOL, and DHC-6 Twin Otter. The 

first three being heavy, transport-type aircraft, and the fourth, a 
light, STOL utility transport. A complete discussion of the aircraft 
equations of motion is included in Appendix A--Equations of Motion 
and expressions for the aerodynamic coefficients are found in 


la 
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Appendix B— Computer Program Documentation. Landing simulations with 
either fixed controls or automatic feedback controls are carried out 
for each aircraft. 


4 



CHAPTER II 


FLIGHT PATH STUDIES — BASIC CONCEPTS 

The flight path simulation studies with fixed controls uses two 
basic non-di mens ionali zed initial altitudes: z = 3.33 (1,000 ft or 

305 m) and z * 1.0 (300 ft or 91 m). The aircraft is trimmed at these 
altitudes for a fixed glide slope angle. The corresponding throttle 
and elevator angle settings are then held constant for the remainder of 
the landing. Different phases of the investigations used different 
glide slope angles in the following manner: 

1. 2.7° for DC-8, B-747: normal ILS glide path. 

2. 6.0° for B-747: steep angle approach. 

3. 7.0° for STOL, DHC-6: normal ILS glide path. 

4. 2.7° for DHC-6: shallow angle approach. 

The introduction of wind shear (gradient of vertical and horizontal 
mean wind components) causes deviation from the glide slope, as well as 
from the expected touchdown points. These deviations were then 
compared to the no-wind condition. From this comparison, the magnitude 
of the deviation was found to be directly related to the associated 
wind shear (negative values for head wind and updraft, positive values 
for tail wind and downdraft). However, some of the values tend toward 
unrealistic effects. These were a direct consequence of the inability 
of the aircraft to negotiate the imposed wind field with the controls 
fixed. Overcoming the unrealistic effects of fixed controls led to the 
development of an automatic control system for the same two-dimensional 
system of equations [6], 
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The first consideration in the creation of an automatic control 
system is to regard every phase of the flight as the accomplishment of 
a set task (i.e., flight along a specified trajectory). Thus, the 
desired path, in this case an ILS beam of 2.7°, is known and the 
departures from it, designated as errors, can be calculated. These 
errors are a result of the imposed wind shear. The control system 
then consists of numerical correction of the errors through thrust, 
elevator angle, and pitch angle commands to maintain a specific 
trajectory. This type of system was tested, in this investigation, 
only for a DC-8 model , 

The control system consists of four phases: hold, capture, 

track, and flare. In this simulation, the DC-8 is initially trimmed 
at 1,000 ft (305 m) for straight and level flight. The first phase 
holds the aircraft's altitude until it intersects the ILS beam. It 
then switches to the glide slope capture mode, consequently activating 
the thrust and elevator controls to capture the beam. As soon as the 
specified glide path is established, the third phase, a glide slope 
tracking mode, becomes effective. The command signals actuate the 
governing control systems such that the DC-8 remains on the glide path. 
This is maintained until switching to the flare mode's initiation 
altitude of 60 ft (18 m) AGL (above ground level). The aircraft 
remains in the flare mode for the duration of the flight (see 
Figure 1). 

For this investigation the scope of the automatic landing 
problem was restricted in two ways: 
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Glide-slope 

capture 



Figure 1. Automatic landing geometry using ILS [6]. 
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1. Aircraft simulation equations are restricted to three 
degrees of freedom by considering the longitudinal axis only, as noted 
in Appendix A--Equations of Motion. 

2. System guidance information is assumed to come from error- 
free sensors and an error-free ILS beam. 

These are reasonable restrictions and may be justified by the 
following: 

1. Wind-shear related accidents from NTSB aircraft accident 
reports indicated a hazard primarily from aircraft being forced below 
minimum safe altitudes, landing short or long, rather than lateral 

variation from the glide slope [6]. 

2. It is beneficial to maintain simplicity of the automatic 
control routines, since the objective of this study is the effects of 
wind shear and not a study of ILS system errors. 
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CHAPTER III 


GENERATION OF THUNDERSTORM AND TURBULENCE MODELS 

The eleven wind fields used in this study are a direct 
measurement of storms passing the WKY-TV/NSSL (National Severe Storms 
Laboratory) 1,500-ft (500 m) meterological tower in Norman, Oklahoma, 
during the months of May to July from 1971 to 1973. The thunderstorm 
wind shear data were recorded in the form of time histories of the 
longitudinal, lateral, and vertical wind vector components measured 
by aerovanes and anemometers [1], Time histories of these data were 
then converted to horizontal spatial coordinates with Taylor's 
hypothesis (i.e., x = Wt), thus forming a two-dimensional wind field 
in a vertical plane. The data were then numerically tabulated in an 
11x41 array. The tabulated wind speed values are coupled with a 
computer look-up routine developed to provide the wind speed values 
and gradients at any position (x,z) called for in the main program. 

A complete discussion of the wind field model is given in [7]. 

Reference [8] notes two potentially dangerous shear situations 
during descent which are inherent to all eleven thunderstorm cases 
studied. The first of these is a tail wind shearing out to a calm or 
head wind component. This results in an increase in the relative air- 
speed and the associated lift causes the aircraft to rise above the 
glide slope (see Figure 2a). To combat this, thrust and/or pitch 
attitude must be continuously reduced. However, an aircraft with 
fixed controls tends to overshoot the touchdown point. 
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(a) Approach in tailwind with shear 



Figure 2. The two main potentially hazardous shear situations [8]. 
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The second hazard situation arises from a head wind shearing 
out to a calm or tail wind component. In this instance, the airspeed 
decreases causing the aircraft to pitch down with an associated loss 
of altitude, and with fixed controls, to land short (see Figure 2b). 

With all thunderstorm gust fronts, strong vertical and hori- 
zontal shears occur [9], Vertical shear is a variation in the wind 
speed components with height z, whereas horizontal shear refers to 
variation with horizontal distance x. In the case of two-dimensional 
wind field vertical updrafts and downdrafts tend to compound the 
effects cited previously. This is explained by the following 
discussion. 

Since the total steady aerodynamic force on an aircraft is 
conventionally decomposed into lift and drag components, the forces 
acting on the aircraft may be resolved, as shown in Figure 3. 

Noting the lift and drag forces always act normal and parallel 
to the relative wind, respectively, the lift and drag in a perturbed 
wind condition is shown in Figure 4. The quantities L Q and D q repre- 
sent the lift and drag forces acting on the airplane during steady 
flight condition, and are normal and parallel to the relative wind — 
taken parallel to the x-axis in this case. Thus, when the aircraft 
encounters an up or down draft, as well as a change in wind speed, the 
relative wind shifts to a new position in a direction opposite to the 

vector sum of V and w , as shown in Figure 4. This perturbed position 
a z 

of the relative wind causes an increase in the angle of attack a p . 

The vector quantities L and D then represent the aerodynamic forces 

acting on the aircraft during the disturbed wind condition [10]. 
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Figure 3. Lift and drag acting on an airplane [10]. 



Figure 4. Variation of lift and drag with change in 
vertical wind component [10]. 


To further illustrate this, consider the simple example of an 
airplane in level flight. The relative airspeed is given by 

- V E , and the lift, C, and drag, 6, are normal and parallel to 
^ a , respectively. Assume |V E | = 100 mph (45 m s' 1 ) and that 
C l/ C D = 2,0j ^ or simplicity* Also assume a 16-mph (7 m s -1 ) tail wind 
shearing to a 4.5-mph (2ms 1 ) head wind, which corresponds to Case 1 
described in the foregoing. Figure 5a shows the vector relationship 
for the initial condition with the forces balanced, and Figure 5b shows 
the latter condition with the resultant force t acting upward and 
backward (forces — subscripted 1). This causes the airplane to go high 
on the glide slope and to increase airspeed. 

If the shear is accompanied by an updraft of 7 mph (3 m s" 1 ) 
or a downdraft of 7 mph (3 m s -1 ) the former case will act approxi- 
mately 4° to the left (subscripted 2), and in the latter case, 4° to 
the right (subscripted 3), as illustrated in the diagram. The 
accompanying updraft creates forces which cause the airplane to be 
higher on the glide slope, whereas the downdraft creates forces which 
cause greater loss of airspeed. 

Recently, McCarthy and Blick [11] performed an analysis for a 
B-727 class airplane with fixed controls encountering a half-sine wave 
tail wind of 23 mph (10 m s -1 ) amplitude at the phugoid frequency along 
a -3° glide slope which experienced an initial decrease in altitude 
from the flight path; this relates to the second type of hazard 
discussed earlier. They found that for any decrease in airspeed V a 

a 

(from an increasing tail wind or decreasing head wind) the lift and 
drag are reduced since there is little change in thrust and 
essentially no change in weight, an unbalanced force forward and down 
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is experienced, causing the airplane to accelerate in this direction 
(see Figure 6). Due to the inertia of the aircraft, the forward 
component of the unbalanced force is not large enough to cause the 
aircraft's speed to increase. This is due to the increased tail wind 
of their study being larger than the inertial airspeed perturbation 
during the first few seconds. However, the inertial airspeed 
eventually overcomes the reduction in the indicated airspeed caused by 
the horizontal gust (tail wind in this case), and the aircraft 
indicated airspeed starts to increase. 

The variations in wind speeds for the thunderstorm model are 
10-sec averaged values (see [7]). The wind field is thus a quasi- 
steady wind field, and to impose gusts of higher frequencies, a model 
of turbulence simulation was used with the random signals being super- 
imposed on the flow field. The random turbulence signal provided by 
the turbulence modeling scheme has a Dryden energy spectrum and 
incorporates the non-Gaussian characteristics of the turbulence [12]. 
Reference [7] provides a complete discussion of the turbulence simu- 
lation technique used in this study. One flight simulation using 
automatic controls employing this turbulence simulation scheme was 
completed for a DC-8 type aircraft. 
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Figure 6. 


Transient condition 
(or increasing tail 


from decreasing head wind 
wind) [11]. 
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CHAPTER IV 

RESULTS AND DISCUSSION 

For the figures appearing in this chapter, numerous computed 
aircraft parameters are plotted as they vary spatially along the 
flight path. However, for some parameters, such as automatic control 
inputs, temporal variations are important, and thus, are shown as 
such. Of principal interest in the fixed control mode is vertical 
position z, horizontal position x, and pitch angle y. Also considered 
is the indicated airspeed (IAS) V . The spatial variations, z and x, 

a 

are non-dimensional i zed with the reference height h = 300 ft (91 m). 

a 

Angle of pitch y is expressed in radians, and V is non-dimensional i zed 

with the initial value of the relative velocity V (i.e., the value 

a 

assigned at trimmed conditions). 

For cases utilizing the automatic control subroutines, the 
magnitudes of the control inputs--thrust Fy and elevator angle 6^-- 
thrust is non-dimensional ized with the thrust at trimmed conditions 
Fj^, which is on the order of 2.87x10^ lbs ( 1 .3x1 0^ kg) for a DC-8. 
Elevator angle is given in degrees, whereas the pitch angle is 
expressed in radians. Vertical and horizontal position of the aircraft 
with automatic controls is given, but only to establish that the air- 
craft maintains the specified trajectory. The control inputs which 
arise from corrections necessary to maintain a given glide slope are 
the variables of interest with an automatic control system, rather than 
deviations from the flight path, as is the case with fixed controls. 
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Table 1 lists the many various flights which have been simulated using 
the computer code documented in Appendix B. 

Fixed Controls 

The nine ill ustrations--Figures 7 through 15 — following Table 1 
constitute the complete set of resulting flight paths for the 120 
simulations with fixed controls. One can see a large variety among the 
particular flight paths traced by the different aircraft, but on the 
large scale a definite recognizable pattern is observed and may be used 
to categorize the flight paths into two overlapping groups. These are 
(1) an oscillatory phugoidal motion (e.g,. Figure 7, DC-8 Cases 1 
through 10; Figure 9, DHC-6 Cases 2, 4, 6, 8, and 9; Figure 12, B-747 
Cases 1 through 10; and Figure 14, augmentor-wing STOL Cases 1, 2, 3, 

5, 6, 9, 10, and 11), and (2) a non-osci llatory motion (e.g., Figure 7, 
DC-8 Case 11; Figure 11, DHC-6 Cases 1 through 11; Figure 13, B-747 
Cases 5, 7, 8, and 11; and Figure 15, augmentor-wing STOL Cases 1 
through 11). As one would expect, there are exceptions to the two 
specified categories, First, in Figure 9, flight through thunderstorm 
Case 10 displays a divergent characteristic different from any other 
flight path. No explanation is offered except that this divergent 
flight path probably is a consequence of the longitudinal dynamics of 
the DHC-6 and the combined head wind and strong updraft conditions 
experienced in the wind field. The second exception can be found in 
Figure 14 with checking the flight paths traced through Cases 7 and 8 
by the augmentor-wing STOL aircraft. Here the aircraft executes a 
loop while operating with fixed controls, possibly due to the air- 
plane's inability to regain stable flight after an out-of-phase wind 





Table 1. 

Flight Simulations 


Ai rcraft 
Type 

Initial 

Conditions 

Special Conditions 

Thunderstorm 

Cases 

A1 ti tude 
(ft Cm]) 

Glide Slope 
Angle 
(deg) 

Ground Automatic 

Effect Controls Turbulence 

DC-8 

1 ,000[305] 

2.7 


1 through 11 


300 [91] 

2.7 


1 through 11 


1 ,000[305] 

2.7 

X 

1,2,9,11 


1 ,000[305] 

2.7 


1,9,11...W Z =0 


1 ,000[305] 

2.7 


1,9,11...W X =0 


1 ,000[305] 

0.0 

X 

1,2,9,11 


1 ,000[305] 

0.0 

X X 

1 

DHC-6 

1 ,000[305] 

2.7 


1 through 11 


1 ,000[305] 

7.0 


1 through 11 


300 [91] 

7.0 


1 through 11 

B-747 

1 ,000[305] 

2.7 


1 through 11 


1 ,000[305] 

2.7 

X 

1 through 11 


1 ,000[305] 

6.0 

X 

1 through 11 

Augmentor- 





wing STOL 

1 ,000[305] 

7.0 

X 

1 through 11 


300 [91] 

7.0 

X 

1 through 11 


X : Indicates special conditions were employed during those runs. 


z/h 





z/h 



Figure 9. Flight paths of DHC-6 landing with fixed controls from 1 ,000-ft (305 m) 
level, glide slope -2.7°. 


Figure 10. Flight paths of DHC-6 landing with fixed controls from 1 ,000-ft (305 m) level, 
glide slope -7.0°. 



Figure 11. Flight paths of DHC-6 landing with fixed controls from 300-ft (91 m) level, glide 
slope -7.0°. 
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Figure 13. Flight paths of B-747 landing with fixed controls from 1 ,000-ft 
(305 m) level, glide slope -6.0°. 





shear. This flight path is, of course, unacceptable and demonstrates 
the STOL aircraft cannot negotiate a thunderstorm with the controls 
fixed. Those flight paths which terminate before landing are a result 
of the aircraft flying outside of the experimental wind field data. 

Since there are two easily distinguishable classes of flight 
paths that have been established, redundancy arising from discussion 
of each individual flight path may be avoided by investigating repre- 
sentative cases of the flight paths. For this study, flight simulation 
data related to thunderstorm Cases 9 and 11 were chosen for detailed 
discussion. However, in some instances for reasons of comparison, 
results from additional cases are supplied. Cases 9 and 11 are 
inherently different in the sense that the aircraft reacts in the 
manner described by the two main hazards in [8] (but not necessarily 
describing the wind conditions causing these effects); therefore, 
their selection for providing representative flight paths. 

Figure 16 illustrates the wind profiles "seen" by a DC-8 while 
landing with fixed controls through thunderstorm Cases 9 and 11. The 
wind speeds, W and W , are non-dimensional i zed with the initial air- 
speed V , A negative value indicates a head wind or updraft, whereas 
a positive value is a tail wind or downdraft for the longitudinal and 
vertical winds, respectively. Of course, the wind profiles shown in 
the figure would be different depending on the point within the flow 
field that the aircraft's flight began. Beginning at different points 
in the wind field also causes a change in the initial trim conditions 
of the aircraft. 

Figure 17 is a comparison of fixed control Cases 9 and 11 with 
the aircraft initially trimmed to follow a 2.7° glide slope. The 
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Figure 17. Comparison of aircraft landing with fixed controls in thunderstorm 
cases 9 and 11 from 1 ,000-ft (305 m) level, glide slope -2.7°. 










aircraft involved are DC-8, B-747, and DHC-6. As to be expected in a 
non-equilibrium condition, none of the flights follow the established 
glide angle. In fact, the deviations shown are quite extreme and tend 
toward intolerable values. For Case 9 the aircrafts’ phugoid oscil- 
lation is wildly excited. The frequency of oscillation for the 
aerodynamically similar DC-8 and B-747 is approximately the same, 
except for a slight phase shift. The oscillation of the slower, 
lighter DHC-6 is less pronounced, but still very much present. The 
non-linear, computer-simulated values and the 1 inear-predicted values 
of the phugoid period and characteristic wavelength are presented in 
Table 2. The predicted values for the phugoid period are given by 
Etkin [13] as T = /l tt V / g. The phugoid wavelength A is then found 

a 

A 

from A = VT, and non-dimensional ized as A = A/h . where h = 300 ft 

a a 

(91 m) . Notice that the simulated values from Figure 17 for a DC-8 
and B-747 correspond closely with the linear-predicted values, but 
this relation does not hold true for the DHC-6. Table 2 also shows 
the phase shift between the DC-8 and B-747, which is dependent on the 
landing speed of the aircraft. 

Tabulated values of the deviation from the expected touchdown 
point presented in Table 3 provide an indication of the severity of 
the wind shear effects due to the proportionality between the wind and 
deviation from the touchdown point, as mentioned previously. Many of 
the flights are forced to land short, a tragic event by any margin, as 
proved by examining wind shear related NTSB aircraft accident reports 
(see [2] through [5]). Overshooting of the touchdown point is not as 
alarming since a go-around can almost always be executed. The vari- 
ability in the range of touchdown points is also indicated by Table 3. 
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Table 2. Phugoid Period and Horizontal Wavelength 



V a 

(ft s' 1 

[m s' 1 ]) 


T(sec) 


X(ft [m]) 



/\ 

A 


Air- 

craft 

Computed 

Pre- 

dicted 

£ 

Computed 

Pre- 

dicted 

£ 

Computed 

Pre- 

dicted 

£ 

DC-8 

230 

[70] 

29.9 

31.7 

1.8 

7,152 

[2,180] 

7,228 

[2,203] 

76 

[23] 

23.84 

24.09 

0.25 

B-747 

217 

[66] 

28.8 

30.0 

1.2 

6,780 

[2,067] 

6,840 

[2,085] 

60 

[18] 

22.6 

22.8 

0.2 

DHC-6 

150 

[46] 

27.1 

20.7 

6.4 

7,890 

[2,405] 

3,333 

[1,016] 

4,556 
[1 ,390] 

26.3 

11.11 

15.19 


GO 

CO 


e: | computed value - linear-predicted value|. 



Table 3. Deviation from Touchdown Point for Condition 1 

(Aircraft Trimmed at 1 ,000-Ft (305 m) Level for 2.7° Glide Slope) 


Air- 





Thunderstorm Case 

Number 





craft 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

n 

DC-8 

OD 

OD 

-1 ,875 
(-572) 

OS 

-3,225 

(-983) 

OD 

-7,425 

(-2,263) 

-8,025 

(-2,446) 

OS 

OD 

-9,150 

(-2,788) 

B-747 

OD 

OD 

-2,100 

(-640) 

OS 

-3,225 

(-983) 

-2,250 

(-686) 

-7,725 

(-2,355) 

-8,475 

(-2,583) 

OS 

-1,800 

(-550) 

-9,450 

(-2,880) 

DHC-6 

-8,550 

(-2,606) 

OD 

-2,850 

(-864) 

+1 ,875 
(+572) 

-5,250 
(-1 ,600) 

-2,550 

(-777) 

-8,700 

(-2,652) 

-7,725 

(-2,355) 

OS 

OS 

-10,350 

(-3,155) 


±n: Distance from touchdown point in feet (meters). 

OS: Severe overshoot--actual value not computed before exhausting data. 

( 

OD: Out of data range. 



To assess the influence of shear in the longitudinal wind 
relative to that in the vertical wind. Figure 18 for Case 9 was 
employed. In this figure, landing with both wind components is com- 
pared to landings with only the individual wind components. Consider 
the line marked “glide slope" as the no-wind condition. This is the 
path the aircraft is initially trimmed to follow. However, after 
introducing a vertical wind component, W z , only, one observes from 
Figure 18 that a slight departure from the glide slope occurs. A hint 
of phugoidal excitation is present also. With only a longitudinal 
wind component, W x , the departure from the glide slope essentially 
coincides with the flight path for the total wind field. This suggests 
that shear in the longitudinal wind is the most significant contributor 
to the excitation of the phugoid mode. This lends support to the con- 
clusion of McCarthy and Blick [14] that the longitudinal wind speed 
wavelength of thunderstorms causes instability in the phugoid mode. 

They have shown that a horizontal gust produces a large peak in 
the aircraft velocity perturbation and a lesser peak in altitude pertur- 
bation at the aircraft phugoid frequency. This means that if a steady 
sinusoidal horizontal gust input on the order of 4 knots (2.06 ms 1 ) 
were encountered, the aircraft, depending on its aerodynamic character- 
istics, would respond with a sinusoidal velocity perturbation of 
approximately 40 knots (20.6 m s” 1 ). At one point in its cycle the 
aircraft would approach a stall speed or go below it during each sine 
wave cycle. These authors found that vertical sinusoidal gusts do not 
affect the aircraft velocity as much as horizontal gusts. They noted 
that three minutes before Eastern flight 66 (Boeing 727) crashed at 

New York's Kennedy Airport on June 24, 1975, a light aircraft 
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(Beechcraft Baron) made a successful landing, although it did experience 
a heavy sink rate and an airspeed drop of 20 knots (10.3 m s"^). Their 
premise is that medium-size jet transport aircraft tend to experience 
larger excursions in velocity and altitude when flying through horizon- 
tal gusts having large spectral components near the phugoid frequency 
than do lighter aircraft. This observation is in complete agreement with 
the results shown in Figure 17. 

On the other hand, Fujita [15] analyzed the same Eastern 66 
accident and attributes the accident to the strong downburst resulting 
from flying through the center of the down draft zone of the thunderstorm's 
cell. Figure 19 shows the flight path of a DC-8 type aircraft landing 
through the wind field associated with the Kennedy accident as tabulated 
for computer application by Keenan [16]. Results for the case of the 
total wind field (flight path A) and for the case where each wind compo- 
nent is individually set equal to zero are shown (flight path B, W x = 0; 
flight path C, W = 0). On top of the figure are shown the wind speeds 
encountered along flight path A. 

The interesting observation is that the downburst of approxi- 
mately W, = 12 m s" 1 at x/h^ = 18.5 applied separately would cause the 
aircraft with fixed controls to crash at approximately x/h g = 30. 

However, when coupled with the increasing headwind, the aircraft manages 
to negotiate the severe down drafts and land although experiencing large 
amplitude oscillations. The horizontal wind shear component alone 
causes less severe flight conditions. 

It appears that the combined effect of both wind shear components 
is important. Had the longitudinal wind speed been shearing out, i.e., 
a decreasing headwind, the aircraft would have landed even shorter. 
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Figure 19 DC-8 type aircraft landing with fixed 
controls in windfield associated with 
JFK Eastern 66 accident and winds 
encountered along flight path A. 
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Therefore, the conclusion of Blick, et al . [17] are not confirmed in 
this study. Further study is required to determine how longitudinal 
and lateral wind shears combine to create hazardous effects. 

Comparing the DC-8 and B-747 flight path angle y and indicated 
airspeed V Q landing through the two representative wind fields 
(Figures 20, 21, and 22) contributes a broadening picture of aircraft 
behavior with fixed controls in thunderstorms. The initial trimmed 
condition of flight path angle (Figure 20) is -0.04712 radian (or -2.7°) 
and is represented by the broken line. Deviation from this line is the 
perturbed flight from the specified glide slope. For Figures 21 and 
22, the indicated airspeed is non-dimensional ized with the initial air- 
craft velocity relative to the air V a , and is plotted against the 

o 

non-dimensional ized horizontal distance x, resulting in a unique speed 

profile for each of the two aircraft considered. Figure 21 indicates 

that the airspeed reaches a low of 118 kts (61 m s - ^) twice for the 

DC-8 model, and two independent lows of 105 kts (54 m s" 1 ) and 100 kts 

(52 m s _1 ) for the B-747. All four low points are attained at a pitch 

angle of 0.0 radian. The level stall speed for a DC-8 and a B-747 is 
-1 -1 

113 kts (57 ms ) and 108 kts (55 m s ), respectively, which indi- 
cates the DC-8 aircraft is operating at an unsafe margin above stall 
and the B-747 would stall. Stall, at the altitudes indicated by 
Figures 7 and 12, pages 20 and 25, respectively, could prove to be 
catastrophic since recovery may not be possible. The flight path angle 
and IAS curves for flight through Case 11 for the two aircraft appears 
more "tame" than for the previous example. Case 9. In fact, the two 
lowest points of Figure 22, representing V 3 p C g = 127 kts (66 m s~^ ) 
and V ag = 120 kts (62 m s - ^), are reasonably above the stall speeds 
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for the respective aircraft. Since these minimum velocities occur at 
an angle other than y = 0, the stall speeds are found by the following 
equations : 

V$dc_q = 147 kts/1.3 cos y or 74 m s“Vl.3 cos y 
and 

V$b _747 = 140 kts/1.3 cos y or 71 m s”Vl-3 cos y [18] . 

Before proceeding to the next figure, a discussion of ground 
effect calculations used in this study is necessary. For the DC-8 
and B-747 the ground effect terms included in the equations of the aero- 
dynamic coefficients are dissimilar in nature, but the effects are 
approximately the same. It was found from flight path studies in 
which 15 comparison flights were conducted (4 for DC-8 and 11 for 
B-747), that inclusion of ground effect had negligible effect on the 
flight path (see Figure 23). However, the ground effect terms for the 
augmentor-wing STOL are an important consideration because they come 
into play at a much higher al titude--close to 200 ft (61 m) AGL. The 
specific ground effect terms used in these computations are given below 
for a DC-8, B-747, and augmentor-wing STOL [19]: 

For DC-8: 

AC. = 0.063(C, ) £a|- , 

L GE L *> GE 

AC d = (-0.02 - 0.332 a')e GE , 

GE 

AC = 0.066(C. )e rf - , 

m GE L oo GE 

where 

e GE = 0.972 e' h/17 , 

h = wheel height , 
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Figure 23. Example of the role of ground effect terms 


C L = free stream value of . 


For B-747 : 

AC, = Kp C (0.240)cos[8.036(a' - 0.00526)] , 

l ge ge 

AC n = K* (2.308 a' 3 - 0.9796 a' 2 - 0.1769 a' - 0.00384) , 
°GE GE 

AC = K G (2.736 a' 2 - 0.621 a' - 0.115) , 
m GE GE 

where 

K* = 1.7034xl0' 6 h 3 - 1.0736xl0" 4 h 2 - 1.4813xl0‘ 2 h + 1.0 , 
GE 

id* = 3.7906X10' 6 h 3 - 4.937xl0~ 4 h 2 + 2.807xl0~ 3 h + 1.0 , 

GE 

noting the ground effect terms are included only during the last 
82.5 ft (25.0 m) of flight. 

The ground effect terms for the ST0L aircraft are more 
important than for conventional aircraft. The complexity of these 
terms is also increased and given by: 


Act 


I = [ C. c.(l - q/qj + C. q/q„][l/ (1 - ^)c ] 

Li in L n J Li in '■'i 


'WB 


GE 


X, 


WB. 


1 WB 


a 


AC n =(C n +C a' ) (q/q - 1 ) - [ C 2 (^)](q/qj 
U GE U o a L WB_ l L 


+ 0 - q/qjc D ex , 
c j 


AC. 


m 


GE 


= ( C + C a' ) (q/q - 1 ) + C ' ( ^ ) C. 

m o m a ” m a C L L WB ge 
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where 


+ C C,(l - q/q ) , 
m r j ^ 

J 


= j 1 + 


WB oo h 2 

^-[ 1 + 16 ] 

8tt h/c c-AR 


- 1/2 


-2 


} • 


and 


2 

#*■ = |2ir AR[ 1 + 16 (r^— ) ]} , 

l L l c-AR > 


C. = thrust coefficient =0.75 (for this study). 
J 


For the DHC-6 ground effect terms were unknown, and consequently, their 
influence was not assessed. 

The influence of glide slope angle on the flight paths through 
the storm is considered next. In Figure 24 the spatial position of the 
aircraft during landing is shown for the augmentor-wing STOL, DHC-6, 
and B-747 in thunderstorm Cases 9 and 11. Notice that two glide 
slopes, 6° and 7°, are indicated in the illustration. The glide slope 
of 7° is the initial trim condition for the STOL and DHC-6 Twin Otter 
and is representative of the usual approach angle for this type of air- 
craft. The 6° glide slope is the initial trim condition for the B-747 
and is included to investigate the effect of approach angle on 
alleviation of the wind shear hazard. This 6° glide slope is not 
unusual for a two-step ILS approach [20]; however, the aircraft usually 
establishes the standard 2.7° approach angle, at a given altitude, for 
the duration of the flight. In this study, however, the B-747 was not 

allowed to make the transition to a reduced approach angle, and thus is 
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(a) Case 9 comparison 



(b) Case 11 comparison 


Figure 24. Flight path comparison of aircraft landing with fixed 
controls from 1 ,000-ft (305 m) level with increased 
approach angle in thunderstorm cases 9 and 11. 



forced to land steep. Landing the aircraft steeply provided infor- 
mation on the effect of increased glide slope angle on conventional 
aircraft. 

The influence of glide slope angle is investigated only for 
Cases 9 and 11 which are indicative of all flights in this phase of the 
study, except for augmentor-wing STOL Cases 7 and 8, which loop at a 
low altitude. As noted previously, this loop characteristic is most 
probably related to the inability of the STOL aircraft to maintain a 
stable condition with fixed controls through these particular thunder- 
storm cases and is not considered a realistic simulation. B-747, 

Case 9, has the most extreme departure from the glide slope, but a 
majority of the B-747 flights follow similar trajectories. The 
DeHavilland Twin Otter's flight path is an appreciably damped version 
of the phugoidal motion experienced for cases where previously the 
glide slope was -2.7°. This is also true, to a reasonable extent, for 
the B-747 (refer to Figure 12, page 25, and Figure 24). This obser- 
vation tends to suggest an inverse relationship between wind shear 
effects and approach angle, as evident from inspection of Table 4. 

Figure 25 illustrates the flight paths for the STOL and DHC-6 
aircraft entering at the 300-ft (91 m) level. Note that the B-747 was 
not run from the 300-ft (91 m) level with a 6.0° glide slope; there- 
fore, it is excluded from the figure. Examination of this figure shows 
that the violent phugoid oscillations present for entry at the 
1,000-ft (305 m) level (see Figure 24) have yielded to a more stable 
flight. Only the augmentor-wing STOL has a remanent oscillatory 
motion. Because of a strong damping effect, brought about by increased 

glide slope angle (-7°) and decreased initial shear encounter altitude 

48 



Table 4. 


Comparison in Touchdown Point Deviations from the 1,000-Ft (305 m) 
Level Using Different Glide Slope Angles 


Thunderstorm Case Number 


Ai rcraft 

1 — 

2 

3 

4 

5 

6 

1 

8 

9 

10 

IT 

Remarks 

Augmentor- 
wing ST0L 

-3,105 

(-946) 

0D 

-2,955 

(-901) 

00 

-2,385 

(-727) 

0D 

LP 

LP 

0D 

-2,265 

(-690) 

-3,405 

(-928) 

Condition 2 

DHC-6 

-270 

+405 

+180 

+480 

-945 

+330 

-1,570 

-1,020 

+330 

+630 

-1,958 

Condition 2 

Twin Otter 

(-82) 

(+123) 

(+55) 

(+146) 

(-288) 

(+101) 

(-479) 

(-311) 

(+101) 

(+192) 

(-597) 



8,280 

(2,525) 

UNK 

2,670 

(814) 

1,395 

(425) 

4,305 

0,312) 

8,370 

(2,550) 

7,170 

(2,045) 

6,705 

UNK 

UNK 

8,302 

(2,560) 

Improvement over 
Condition 1; 
n:distance closer 
to T.D. point 

B-747 

-690 

(-210) 

+1,290 

(+393) 

+840 

(+256) 

+1 ,770 
(+540) 

-810 

(-247) 

+450 

(+137) 

-1,650 

(-503) 

-660 

(-202) 

+690 

(+210) 

+1 ,470 
(+448) 

-2,190 

(-668) 

Condition 3 

WO 

UNK 

UNK 

1 ,260 
(384) 

UNK 

2,415 

(736) 

1 ,800 
(550) 

6,075 

(1,852) 

7,815 

(2,382) 

UNK 

330 

(100) 

7,260 

(2,213) 

Improvement over 
Condition 1; 


n:distance closer 
to T.D. point 


±n: Deviation from touchdown point in feet (meters). LP: Looped. 

OD: Out of data range. UNK: Difference not calculable, previous value unknown. 
Condition 1: Aircraft trimmed at 1,000-ft (305 m) level for -2.7° glide slope. 
Condition 2: Aircraft trimmed at 1,000-ft (305 m) level for -7.0° glide slope. 
Condition 3: Aircraft trimmed at 1,000-ft (305 m) level for -6.0° glide slope. 



300 ft (91 m) 


DHC-6 


Glide 
S 
/ 



Case 11 comparison 


Flight path comparison of STOL and DHC-6 landing with 
fixed controls from 300-ft (91 m) level, glide slope 
-7.0°, in thunderstorm cases 9 and 11. 




(z = 300 ft = 91 m) , many of the flights land within close proximity 
to the desired touchdown point, as shown in Table 5. 

Table 5 shows the influence of height at which the aircraft 
enters the wind field on the departure from expected touchdown. It is 
shown that less deviation in touchdown occurs for entry at a lower 
altitude than for entry at a .higher altitude. Of course, this effect 
is strongly dependent on the characteristics of the storm and the 
particular point of entry. These deviations are felt to be repre- 
sentative since eleven arbitrary thunderstorms are considered. 

However, considerably more simulations starting at different positions 
in the flow field are necessary before any positive conclusions can be 
drawn . 

Automatic Controls 

Figure 26 illustrates several of the test simulations using 
fixed controls, automatic controls, and flight with turbulence with 
airplane characteristics of a DC-8. Since the fixed control flight 
paths have already been reviewed in previous text and illustrations, 
this discussion is devoted to the automatic control landing cases. 

From the figure it is seen that Cases 1, 2, and 9 remain on the desired 
ILS glide slope until flare initiation altitude, at which time the DC-8 
prepares to land. On the other hand. Case 11 assumes an altitude 
slightly below the glide slope--a condition arising from the combi- 
nation of sustained head wind and severe vertical downdraft 
sufficiently strong that the automatic landing system could not over- 
come it. However, the DC-8 does flare to land "on target." For the 
simulation having turbulence effects included (Case 1), the resulting 
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Table 5. Comparison in Horizontal Deviation from Touchdown Points from 
Different Altitudes with Different Glide Slope Angles 


Thunderstorm Case Number 

Aircraft 1 2 3 4 5 6 7 8 9 10 ll Remarks 


Augnentor- 
wing ST0L 

-417 

(-127) 

-132 

(-40) 

-837 

(-255) 

-12 

(-4) 

-942 

(-287) 

+258 

(+79) 


2,688 

(819) 

UNK 

2,118 

(645) 

UNK 

1 ,443 
(440) 

UNK 

DHC-6 

Twin Otter 

-402 

(-123) 

+1,533 

(+467) 

+78 

(+24) 

-72 

(-22) 

+33 

(+10) 

-72 

(-22) 


-132 

(-40) 

-1,128 

(-344) 

102 

(31) 

408 

(124) 

912 

(278) 

258 

(79) 

cn 

ro 

8,148 

(2,485) 

UNK 

2,772 

(845) 

1,803 

(550) 

5,217 

(1,590) 

2,478 

(755) 


-282 

(-86) 

-132 

(-40) 

-567 

(-173) 

-792 

(-241) 

-102 

(-31) 

Condition 4 

UNK 

UNK 

UNK 

1,473 

(450) 

2,943 

(897) 

Improvement over 
Condition 2; 
n:distance closer 
to T.D. point 

-162 

(-49) 

-162 

(-49) 

-36 

(-11) 

-72 

(-22) 

-342 

(-104) 

Condition 4 

1 ,368 
(417) 

858 

(262) 

294 

(90) 

558 

(170) 

1,616 

(493) 

Improvement over 
Condition 2; 
nrdistance closer 
to T.D. point; 
minus values indi- 
cate increased 
landing distances 

8,538 

(2,600) 

7,563 

(2,305) 

UNK 

UNK 

10,008 

(3,050) 

Improvement over 
Condition 1 


±n: Deviation from touchdown point in feet (meters). 

UNK: Difference not calculable, previous value unknown. 

Condition 1: Aircraft trimmed at 1 ,000-ft (305 m) level for -2.7° glide slope. 

Condition 2: Aircraft trimmed at 1 ,000-ft (305 m) level for -7.0° glide slope. 

Condition 4: Aircraft trimmed at 300-ft (91 m) level for 7.0° glide slope. 



Figure 26. Flight path comparison of DC-8 landing with (1) fixed controls, (2) 

automatic controls, and (3) automatic controls with turbulence included 
in several different thunderstorm cases. 




aircraft motion exhibits a mild oscillation, due to the variability of 
the wind field, and also immediately departs from the desired 
trajectory. It should be noted, however, that each flight path will be 
different when turbulence is superimposed because of the random nature 
of the flow field. At approximately x = 41.5 the DC-8 flight simu- 
lation terminates due to flight outside the wind data provided for 
thunderstorm Case 1 . 

Figures 27 and 28 illustrate typical wind profiles along a 
-2.7° glide slope for thunderstorm Cases 9 and 11, respectively. These 
wind profiles, associated with the DC-8 having automatic controls, 
appear considerably different from the same thunderstorm cases profiled 
for the DC-8 in Figure 16, page 30, since the path traced by the air- 
craft is different. Notice that both Cases 9 and 11 show similar 

longitudinal head wind components W , but the vertical wind components 

X 

A 

W z for the two cases differ (particularly below z < 2.0). Again, the 
wind speeds have been non-dimensional ized with the initial airspeed V . 

Next, consider Figure 29 which shows the thrust command signal 
FT as a function of time T for the DC-8 with automatic controls. 

FT has been non-dimensional ized by the equivalent thrust FTD for which 
Case 9, FTD = 2,65 x 10 4 lbs (1.29 x 10 4 kg), and for Case 11, 

FTD = 1.87 x 10 4 lbs (8.43 x 10^ kg). From the literature it is 
quite plain that frequent and abrupt changes in thrust are necessary 
to maintain the specified glide slope. In fact, the time required to 
increase the thrust needed to maintain the desired position is often 
greater than the actual rate at which thrust can be increased for the 
DC-8 engines when encountering shear conditions. The rate at which 


54 






CJ1 



f 

.0 L 


10 


20 


30 40 50 60 70 

T (sec) 


Figure 29. Rate of change of thrust required of DC-8 landing with automatic control 
system in thunderstorm cases 9 and 11. 


thrust can be increased for the Pratt-Whitney JT-3D engines on a DC-8 
is given approximately by 7.92 x 10 3 lbs s ^ (3.6 x 10 3 kg s ^ ) [21]. 
Changes in thrust beyond the capability of the engine are not possible, 
thus precluding the possibility of the aircraft's remaining on the ILS 
glide slope. Because this maximum thrust exceedance occurs on certain 
excursions, a thrust rate limiter was added to the automatic control 
system. The approach paths followed with thrust rate limiters being 
utilized are compared with those having no thrust rate limitation in 
Figure 30. The flight path of the DC-8 type aircraft without any limit 
on the maximum thrust, departure from the desired trajectory is approx- 
imately 29.5 ft (9 m), whereas with the limiter, the maximum departure 

is 65.6 ft (20 m). The thrust rate limited aircraft does not follow ILS 

beam as closely in the earlier part of the approach, but both simulations 

do intercept the beam at approximately x = 10,500 ft (3200 m), and 
track precisely along the beam flaring on target. 

The other important command signal, elevator angle <5^, is 
represented graphically in Figures 3] and 32 as a function of time. 

Here, reasonable aircraft response to rate change of the elevator angle 
to maintain the glide slope is demonstrated by the DC-8 within thunder- 
storm Cases 9 (Figure 31) and 11 (Figure 32). Observation shows the 
maximum d6^/dt occurring in either of the two illustrations is 


dS, 


dt 


= 3 .3° s 


-1 


max 


in Figure 31. This is well within the range of response for a DC-8, 
since aerodynamic response times are quite small — on the order of 
10~ 3 sec [22]. 
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Figure 30 Comparison of flight paths for DC-8 type aircraft with and without 
thrust rate limiter. 
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Figure 31. Rate of change of elevator angle required of DC-8 landing wi 
control system in thunderstorm case 9. 
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The previous discussion concerning automatic control systems is 
actually an ‘'idealized" system. In this system, the control loop 
assumes that the variables a 1 , a 1 , 6, and q can be monitored during the 
approach and fed back into the landing system to calculate the variable 
gains continuously along the glide slope. The idealized control loop 
also assumed that the ground speed components z/v and x/v are available 
as feedback inputs to the sevro-mechanisms for the thrust and elevator 
control . 

An alternate or simplified automatic control loop assumes that 
only the relative airspeed, V . is monitored during the approach and is 

Q 

available for computing the variable gains. Additionally, this control 
loop does not allow for z and x as feedback input, but rather expresses 
z/V as -sin y and x/V as cos y. The value of y was set equal to zero 
during the hold and flare modes and then equal to the glide path angle 
during the capture and tracking modes. 

The results of this more "conventional" landing system as opposed 
to the "idealized" control system may be found in Figure 33. The figure 
shows a DHC-6 type aircraft landing through thunderstorm Case 9 wind 
field. Neither simulation establishes the designated trajectory along 
the ILS beam. This is probably due to inadequate specification of 
capture control parameters for the DHC-6 type aircraft. However, the 
tracking control does establish a well defined trajectory which includes 
a successful landing. The results shown in Figure 33 are indicative of 
the control which can be achieved through rather severe thunderstorm 
wind shear with appropriate control systems. 
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a 


Figure 33 Comparison of flight paths with idealized and more conventional control 
systems for DHC-6 type aircraft. 



CHAPTER V 


CONCLUSION 

From the results presented in this investigation it is evident 
that thunderstorm outflows are a hazard to aircraft operations in and 
around the airport terminal area. The danger primarily exists here, 
rather than during cruise flight, because of the aircraft's particular 
configuration (namely, little speed margin above stall, high drag, and 
for some airplanes, limited excess power) and because of its close 
proximity to the ground. 

It has been shown that serious departures from the glide slope 
occur during landing simulations of aircraft having fixed controls 
(i.e., no pilot input) through thunderstorm gust fronts. Corre- 
spondingly, control inputs are large and response times are small for 
aircraft using an automatic control system in the same thunderstorm 
cases. 

The results of this study show that the phugoid oscillations of 

aircraft with fixed controls landing through typical thunderstorm gust 

fronts are highly amplified. This is particularly true for the larger 

transport type aircraft. The amplitude of the oscillations tend to be 

reduced for lighter type aircraft typical of a DHC-6 Twin Otter. This 

is partly due to the characteristics of the aircraft and partly due to 

the steeper landing paths followed. The larger transport type aircraft 

approaching at steeper glide paths than those conventionally used in 

aviation operations have somewhat damped phugoid oscillations, but 

still experience large excursions from the desired landing path. The 
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strong influence of horizontal gradients in thunderstorm associated 
wind fields on the phugoidal oscillation support the Blick, et al . [ 15 ] 
hypothesis that accidents associated with commercial aircraft landing 
through thunderstorm gust fronts may result from the horizontal wind 
shear, but do not rule out that severe downburst can be equally respon- 
sible for accidents. Investigation of flight through the thunderstorm 
wind fields established by Fujita [16] indicate that the downburst can 
cause accidents. On the other hand, the downburst combined with the 
longitudinal wind shear, for the one case studied results in the air- 
craft negotiating the wind field. However, it does experience severe 
oscillations . 

Automatic control systems using variable gains can almost com- 
pletely eliminate the severe perturbations from the flight path for the 
11 thunderstorm models considered in this study. This does not imply 
that automatic control systems can be utilized to land aircraft through 
thunderstorms in all situations. The thunderstorm models utilized in 
this study are obviously not all inconclusive and represent only the 
gust front portion. The automatic control systems have not been applied 
to the extreme downburst winds that have been reported to occur in the 
center of the thunderstorm cell [16]. Moreover, the computer simulation 
treats only two-dimensional effects and therefore excludes additional 
control inputs required to stabilize roll and yaw motions. 

The fact that the automatic control systems do however appre- 
ciably eliminate flight path excursions tend to support the arguments 
that accidents in thunderstorms are not a result of aircraft limitations 
but often are precipitated during transition from automatic to manual 
control [4]. 


65 



REFERENCE LIST 


66 



REFERENCE LIST 


1. Goff, R. Craig. "Thunderstorm Outflow Kinematics and Dynamics," 

NOAA Technical Memorandum ERL NSSL-75, National Severe Storms 
Laboratory, Norman, Oklahoma, December, 1975. 

2. National Transportation Safety Board. "Eastern Airlines, Inc., 

Boeing 727-225, John F. Kennedy International Airport, 

Jamaica, New York, June 24, 1975," Aircraft Accident Report 
No. NTSB-AAR-76-8, National Transportation Safety Board, 
Washington, D. C., March 12, 1976. 

3. National Transportation Safety Board. "Continental Airlines, 

Inc., Boeing 727-224, N88777, Stapleton International Airport, 
Denver, Colorado, August 7, 1975," Aircraft Accident Report 
No. NTSB-AAR-76-14, National Transportation Safety Board, 
Washington, D. C., May 5, 1976. 

4. National Transportation Safety Board. "Iberia Lineas Aereas de 

Espana (Iberian Airlines), McDonnell Douglas DC-10-30, EC 
CBN, Logan International Airport, Boston, Massachusetts, 
December 17, 1973," Aircraft Accident Report No. NTSB-AAR- 
74-14, National Transportation Safety Board, Washington, 

D. C. , November 8, 1974. 

5. National Transportation Safety Board. "Delta Airlines, Inc., 

Douglas DC-9-32, N3323L, Chattanooga Municipal Airport, 
Chattanooga, Tennessee, November 27, 1973," Aircraft 
Accident Report No. NTSB-AAR-74-1 3, National Transportation 
Safety Board, Washington, D. C., November 8, 1974. 

6. Reddy, K. R. "Investigation of Aircraft Landing in Variable Wind 

Fields." Unpublished Master's thesis. The University of 
Tennessee, Knoxville, 1977. 

7. Frost, W., D. W. Camp, and S. T. Wang, "Wind Shear Modeling for 

Aircraft Hazard Definition." Report prepared under Inter- 
agency Agreement No. D0T-FA76-WAI-620, supported by Federal 
Aviation Administration, Wind Shear Office, by FWG Associates, 
Inc., Tullahoma, Tennessee, December, 1977. 

8. Low Level Wind Shear . Advisory Circular Account No. 00-50. 

Washington, D. C. : Department of Transportation, Federal 

Aviation Administration, 1976. 

9. Frost, W., and D. W. Camp. "Wind Shear Modeling for Aircraft 

Hazard Definition." Interim report No. FAA-RD-77-36 prepared 
for U. S. Department of Transportation, Federal Aviation 
Administration, Systems Research and Development Service, 
by FWG Associates, Inc., Tullahoma, Tennessee, March, 1977. 

67 



10. 


McRuer, D. , I. Ashkenas, and D. Graham. Aircraft Dynamics and 
Automatic Control . Princeton, New Jersey: Princeton 

University Press, 1973. 

11. McCarthy, J., and E. F. Blick. "Effects of Wind Turbulence and 

Shear on Landing Performance of Jet Transports." Paper 
presented at the AIAA 16th Aerospace Sciences Meeting, 
Huntsville, Alabama, January, 1978. 

12. Reeves, P. M., R. G. Joppa, and V. M. Ganzer. "A Non-Gaussian 

Model of Continuous Atmospheric Turbulence for Use in Air- 
craft Design," National Aeronautics and Space Administration 
CR-2639, Ames Research Center, Sunnyvale, California, 

January, 1976. 

13. Etkin, B. Dynamics of Atmospheric Flight . New York: John Wiley 

and Sons, Inc., 1972. 

14. McCarthy, J., and E. F. Blick. "Aircraft Response to Boundary 

Layer Turbulence and Wind Shear Associated with Cold-Air- 
Outflow from a Severe Thunderstorm. " Paper presented at the 
7th Conference on Aerospace and Aeronautical Meteorology and 
Symposium on Remote Sensing from Satellites, American Meteor- 
ological Society, Melbourne, Florida, November 16-19, 1976. 

15. Blick, E. F., et al . "Effect of Wind Turbulence and Shear on 

Landing Performance of Jet Transports . " Paper presented at 
the AIAA Conference, Huntsville, Alabama, 1978. 

16. Fujita, T., and F. Caracena. "An Analysis of Three Weather Related 

Aircraft Accidents." Research conducted at the Department of 
the Geophysical Sciences, The University of Chicago, Illinois, 
1977. 

17. Keenan, M. G. Personal communications, October 1977. 

18. Bliss, Jack. Private communication. Flying Tiger Airlines, 

Los Angeles, California, November, 1977. 

19. Luers, J. K. , and J. B. Reeves. "Effect of Shear on Aircraft 

Landing," National Aeronautics and Space Administration 
CR-2287, George C, Marshall Space Flight Center, Huntsville, 
Alabama, July, 1973. 

20. Horonjeff, R. Planning and Design of Airports . Second edition. 

St. Louis: McGraw-Hill, Inc., 1975. 

21. Smith, V. K. Private communication. The University of Tennessee 

Space Institute, Tullahoma, Tennessee, October, 1977. 


68 



22. Norman, W. S. Unpublished classroom notes. The University of 

Tennessee Space Institute, Tullahoma, Tennessee, November, 
1977. 

23. Frost, Walter. "Flight in Variable Mean Winds." Report in 

preparation under NASA Contract No. NASQ-29584 by The 
University of Tennessee Space Institute, Tullahoma, Tennessee, 
1977. 

24. Ralston, A., and H. S. Wilf. Mathematical Methods for Digital 

Computers . New York: John Wiley and Sons, Inc., 1960. 


69 


appendices 


70 



APPENDIX A 


EQUATIONS OF MOTION' 


Equations of Motion 

The two-dimensional model for aircraft motion presented in this 
appendix follows the general form developed by Frost [23]. It accounts 
for both vertical and horizontal mean wind components having both time 
and spatial variations. 

The aircraft trajectory model employed in this study was 
derived based on the following assumptions: 

1. The earth is flat and non-rotating. 

-2 

2. The acceleration of gravity, g, is constant, 32.2 ft sec 
(9.8 m s' 2 ). 

3. Air density is constant, 0.002378 slug ft 2 (12 kg m 2 ). 

4. The airframe is a rigid body. 

5. The aircraft is constrained to motion in the vertical plane. 

6. The aircraft has a symmetry plane (the x-z plane). 

7. The mass of the aircraft is constant. 

8. Initial flight conditions are for steady-state flight. 

Figure 34 illustrates the forces acting on the aircraft. These 

include: 

F-j. = thrust of the engines, 

L = lift, 

-y 

D = drag , 
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FRL 



Figure 34. Relationship between the various forces acting on an 
aircraft [23], 
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W = wind vector. 


mg = gravitational force. 

The figure shows the orientation of the forces with respect to 

->■ 

the velocity relative to the earth (V), the velocity relative to the 
air mass (VK and the fuselage reference line ( FRL ) of the aircraft. 

a 

The x-axis is parallel to the surface of the earth and the z-axis is 
perpendicular to the surface of the earth (positive downward). 

From a direct force balance along the direction V and along 
the direction perpendicular to V, respectively, it follows from 
Figure 34 that 

mV = - L sin 5 - D cos S - mg sin y + Fy cos (6-j. + a) (1) 

and 

mV y = L cos <5 - D sin 6 - mg cos y + F-j. sin (Sy + a) . (2) 

The aerodynamic forces and the thrust from the engines exert a pitching 
moment on the aircraft. The equation describing the momentum balance 
about y is 



F T L T . JML 

I *1 * 

yY *YY 


(3) 


where the dot refers to the derivative with respect to time, and 

g = magnitude of the acceleration of gravity, 

V = magnitude of the velocity relative to the earth, 

y = angle between V and x-axis (the flight path angle), 

F t = magnitude of the thrust vector, 
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m = aircraft mass, 

6y = angle between the thrust vector and the fuselage 
reference line (FRL), 

a = angle between V and the FRL, 

(5 = angle between V* and V, 

q = time derivative of the pitching rate (q), 

lj = effective moment arm of the thrust vector, 

M = pitching moment, 

I yy = moment of inertia about the symmetry plane of the 
ai rcraft. 


By considering a different coordinate system in which the 

x-axis is along the vector V. called "wind" frame of reference by 

Etkin [13], similar force equations can be developed by summing up the 

forces parallel and perpendicular to V . These are 

a 


m ( V + W ) + m q W 
a x w w z w 


= F- 


- D - mg sin y' 


(4) 


w 


and 


m W z - m q w ( V a 
w 


+ W v ) = F t - L + mg cos y 


(5) 


w 


"w 


It is convenient to express these in terms of the wind com- 
ponents relative to an earth-fixed coordinate system, since most wind 
correlations from the meteorological literature are expressed in such 
coordi nates . 
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( 6 ) 


W = W cos y* - W sin y 1 , 
x w x z 


W = W sin y' + W cos y' • 
z w 


(7) 


Taking the time derivative of W , 


w 


W v = W v cos y‘ - sin y' - W v sin y 1 - w z cos y' ^ 


x w x 


Then, since 


dY* 

q w = It ’ 


W = W cos y 1 - W sin y‘ - W q 
*w * w 


(9) 


Also, since 


and 


F^ = Fj cos (6^ + a' ) 
x w 


F-j- = Fy sin(6y + a' ) , 
Z w 


Eq. (4) becomes 


m V a = F t cos(6 t + a 1 ) - D - mg sin y 


- m(W cos y' - W sin y' ) • 

A L 


( 10 ) 


From Eq. (7), taking the time derivative of W , 


w 


W z = W x sin y ' + W z cos Y 1 + W x cos y‘ — W z sin y' ^ 
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and Eq. (5) becomes 


- m q w = - Fj sin(6y + a' ) - L + mg cos y* 

- m(W sin y' + W cos y 1 ) . (12) 

X z 

The moment equation remains the same as Eq. (3). The governing 
force equations in “wind" frame of reference are thus 

m V g = Fy cos (6-j- + a* ) - D - mg si n y 1 

- m(W cos y' - W sin y ' ) , (13) 

X z 

m y' = Fy sin(6y + a‘) + L - mg cos y' 

• * 

+ m(W sin y 1 + W cos y') , (14) 

X 4 


and 


• FfLf 

q = — ■ — - + 
M I YY 


YY 


(15) 


where 


= horizontal wind speed, 

W z = vertical wind speed, 

a' = (angle of attack) angle between V and the FRL. 

a 


Incorporation of Wind Shear 

The wind is seen to enter the equations in the form of a 
gradient or wind shear li x and W z . The expanded form of these 
equations is: 
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J 


!\ + !^cix + 5ldz 

x 9t 9x dt 9z dt 


W„ = 


or 


9W 


W.. = 


9W 


9W 


x at 


r + v [ cos y - sin y -t~] , 


9z 


( 16 ) 


and, similarly. 


W_ = 


z 9t 


9W 9W 9W 

z + V [ cos y - sin y 


9x 


9z 


07) 


Thus, both spatial variations and temporal variations in atmospheric 
motion influence the equations in the wind coordinate system. 

Generally, care is needed in evaluating 9W /9z and 9W /9z, 
since the wind speed is normally expressed in terms of altitude 
measured upward from the surface of the earth, whereas in aerodynamic 
coordinates, z is measured downward. 

Additional kinematic relationships necessary to solve for the 
aircraft motion are as follows: 

The relative velocity as a function of inertial and wind 
velocity is 

o o 1/2 

V a =[(X - W x r + (Z - w z ) 2 ] , (18) 

and, in turn, 

2 

V = W cos y - W sin y + E (W_ sin y - W cos y) 

X Z Z X 

+ v 2 - (W 2 + W 2 ) ] 1/2 . (19) 

a x z 
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-y -y 

The angle between V and is given by 

W sin y + W cos y 

sin 6 = — — - . (20) 

v a 

Other angular relationships are: 

a'=0-Y-d=e-Y l (21) 

and 

a = 6 - y • (22) 

The derivative of a 1 is 

• * • • . 

a 1 = e - Y = q - Y 9 

where y‘ is given by Eq. (14), hence. 



q - 


F-p sin(6y + a' ) 
m V_ 

a 


L 

m 

a 


+ 9 C P S - T ., + J-[W^ sin y' + W z cos y 1 ] • (23) 

v a v a 

Also required for solution of the preceding equations are the 
aerodynamic coefficients: 

C. = C|^ ( ct 1 , 6^, V, q, ol ) , 

Cq = Cp)(a , dp-, V, q, ot , C^) » 

C m = c m ( a ' 9 V, q, a' ) , (24) 

where 

6 E = elevator deflection angle. 
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As is indicated above, the aerodynamic coefficients are 
functions of a number of variables. The expressions for C L , C D , and 
C m , along with the stability derivative data and aircraft physical data 
are given in the following sections. 

The equations of motion discussed in this appendix can be 
solved for the flight of an aircraft flying through spatially and 
temporally varying two-dimensional wind fields. 
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APPENDIX B 


COMPUTER PROGRAM DOCUMENTATION 


Program Description 

This part of the appendix describes a computer program, written 
in FORTRAN IV, designed to solve the equations of motion for aircraft 
flight through spatially and temporally varying two-dimensional wind 
fields developed in Appendix A. The aircraft may assume a fixed 
control mode (i.e., neglects pilot input) or use the automatic control 
subroutines which provide adjustments to the control surfaces (see 
Figures 35 and 36), In addition, an option to include a non-Gaussian 
turbulent field with the modeled thunderstorms for the automatic 
control mode is described schematically in Figure 37 and detailed 
in [7], 

Main Program 

The main calling program first reads the aircraft physical data 
and stability derivatives from the user-supplied data deck. Next, the 
dimensionless groups are calculated. The main program then calls sub- 
routine I N IT which supplies the initial conditions for the equations of 
motion. Subroutine EQUIL, which computes initial values for thrust, 
elevator angle, and angle of attack from the equilibrium conditions, is 
called secondly. Now with initial conditions, error weights , and so 
forth, for the Runge-Kutta scheme defined, subroutine RKGS is called 
and does not return to the main program until the aircraft has landed. 
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Figure 35. Main calling program flow diagram for aircraft 
with fixed controls. 
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Figure 36. Main calling program flow diagram for aircraft with 
automatic controls. 
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Figure 37. Main calling program flow diagram for aircraft with 
automatic controls, including turbulence generation 
subroutines. 
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Nomenclature 


Main program input data: 


NCR 


Read unit number 

NWR 


Write unit number 

PI 

7T 

PI = 3.1415927 

RHO 

P 

RHO = 0.075 

G 

g 

Acceleration due to gravity 

TEST 


TEST = 0 

NP 


NP = 0 

NCO 


NCO = 20 

QV 

qV 

QV = p/2 

PLANE 


Airplane name 

AMASS 

m 

Aircraft weight (lb) 

SA 

s 

p 

Aircraft wing area (ft ) 

CORD 

c 

Wing mean aerodynamic chord (ft) 

VAIR 

V a 

Airspeed (ft s”^ ) 

ALT 

L t 

Effective moment arm of thrust vector 

HA 

h 

a 

Reference altitude, HA = 300 ft (91 m) 

YYI 

Iyy 

Moment of inertia about the symmetry plane of the 
ai rcraft 

DELT 

6 I 

Angle between the thrust vector and the fuselage 
reference line 

CLO 

S, 

Lift parameter 

CLA 


Lift parameter 

CLDEL 

V 

Lift parameter 


84 



I 


CLQ 

% 

Lift parameter 

CLADOT 

C L- 

a 

Lift parameter 

CDO 

C °o 

Drag parameter 

CDA 

C D 

a 

Drag parameter 

CDA2 

C D 2 
a 

Drag parameter 

CMO 

c % 

Pitching moment parameter 

CMA 

C 

m 

a 

Pitching moment parameter 

CMDEL 

C m 
rn r. 

S E 

Pitching moment parameter 

CMQ 

\ 

Pitching moment parameter 

CMADOT 

C m- 

a 

Pitching moment parameter 


Main | 

Drogram dimensionless groups: 

DN1 


QV*SA*HA 

AMASS 

DN2 


G*HA 

(VAIR) 2 

DN4 


CORD 

2HA 

DN5 


QV*SA*C0RD*HA 2 

G*YYI 

DN6P 


HA*G 

AMASS *V AIR 2 

DN7P 


ALT*HA 2 

(vair) 2 *yyi 

PRMT 


Input-output vector for subroutine RKGS 

PRMT ( 1 ) 


Lower bound of the interval PRMT(1)=0.0 
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PRMT ( 2 ) 
PRMT(3) 

PRMT (4) 

DERY 

NDIM 

Listing 


c 

c 

1 

33 


Upper bound of the interval PRMT(2)=100 
Initial increment of the independent variable 
PRMT(3)=0.01 

Upper error bound PRMT(4)=1.0 

Input vector of error weights for subroutine RKGS 

Number of equations to be integrated NDIM=6 

of Main Program 


DIMENSION PLANE ( 20 ) 

DIMENSION A(3,3),BC3) 

DIMENSION Y(6) , DERY ( 6) ,AUX(8,6),PRMI(5) 

DIMENSION DXX(2,23 ,DXZC2,23 , DZX(2,2) , DZZ(2,23 

CDMMON CL,CLO,CLA,CLDEL,CLQ,CLADOT,CLGE 

COMMON CD,CDO,CDA,CDA2, CD3E 

COMMON CM,CMO,CMA, CMDEL,CMQ ,CMADOT, CMGE 

COMMON QV, AMASS,G, Ply S A, CORD , ALT, HA, YYI, VAIR 

CDMMON DN1 ,DN2,DN3,DN4,DN5,DN6,DN7 ,DN6P,DN7P 

COMMON FT, DELI, DEL, DELE, VA, GAMP 

COMMON TEST ,NS,NP,NCQ 

COMMON/UN KOWN/Y 

COMMON /VARIBL/ DERY 

C3MM0N/INI/V0 

COMMON/ WINDS /WX , WZ, *1X1, XZT , WXX , WXZ , WZX , NZZ , WXDOT , WZDOT , ZO , USTAR 
CDMMON/TT/AXA C41 , 1 1 , 2) , DX , DZ 
EXTERNAL FCT , OUTP 

NCR = 5 
NP = 0 
NWR = 6 
NCD = 20 

PI = 3. 141S9265 
KHO=0. 07S 
G=32 . 2 
QV = 0. 5*RHO 

RUNGE-KUTTA PARAMETER SET UP 
NDIM=6 

DO 1 1=1, NDIM 
DERY(1)=0.001 
PRMT (13=0. 

PRMT(2)=100. 

PRMT(33=0.01 
PRMT ( 43 =1 . 0 
TES T=PRMT ( 1 ) 

READ (NCR, 333 PLANE 
FORMAT ( 20 A4 3 
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WRIXE(NWR,41)PLANE 

READ (NCR, 34) AM ASS , S A , CDRD , VAIR , ALT , H A , Y Y1 , DELT 
41 F0RMAT(6X,20A4) 

34 FORMAT (8F10.5) 

READ (NCR, 35) CLD , CL A , CLDEL , CLG , CLADOT 

35 F0RMAT(5F10. 5) 

READ ( NCR , 36 ) CDO,CDA,CDA2 

36 FORMAT ( 3F1 0 . 5 ) 

READ (NCR, 37) CMO , CM A , CM DEL , CMQ , CM ADD! 

37 FORMAT< 5F1 0 . 5 ) 

RRITE(NWR,43) 

43 FORMAT ( 80 ( 1 S 1 ) ) 

WRITE(NWR,44) AMASS, SA, CORD, VAIR, ALT, HA , YYI , DELT 

44 FORMAT ( 6X , ’AMASS = • , FI 0 . 2 , / , 6X , ' h' I NG AREA = » , F 1 0 . 4 , / , 6X , 

1 ' CORD = » ,F10.4,/,6X, ' VAIR = ', F 1 0 . 4 ,/, 6X ALT ( MOMENT ARM) = ', 

2F 8.4,/,6X,«HA = * , FI 0 . 4 , / , 6X , • Y Yl= ' , El 3 . 5 , / , 6X , ' DELT =',F10.5) 
DELI=DELT*PI/180. 

MRlfE(NWR,42 ) 

42 FORMATdX, ' LIFT COEFFICIENT ' ) 

WRITE ( N WK , 4 5 ) CLD , CL A , CLDEL , CLQ , CLADOT 

45 FORMATdX, 1 CLO = 1 , FI 0 . 5 , ' CLA= • , F 1 0 . 5 , « CLDEL= ' , FI 0 . 5 , 

1 ' CLO= ' , FI 0.5, • CLADOT= • , FI 0.5) 

WSlTE(NWR,4t>) 

46 FORMATdX, ' DRAG C0EF1 C I EN IS • ) 

WRIIE(NWR,47) CD0,CDA,CDA2 

47 FORMAT ( 2X , » CDO= ' , FI 0 . 5 , ' CDA= « , FI 0 . 5 , ' CDA2= « , FI 0 . 5 ) 

WRITE(NWR,48) 

48 FORMATdX, ' MOMENT C0EF1 Cl ENTS ' ) 

WRITE(NWR,49) CMO, C^A, CMDEL,CMQ,CMADGT 

49 FORMAT (2X, ' CMO= ' , F 1 0 . 5 , ' CM A = ' , F 1 0 . 5 , 1 CMD£L= • ,F10.5, 'CMQ=' ,F10.5, 
1 ' CM ADOT= 1 ,F10.5,/) 

C CALCULATE DIMENSIONLESS GROUPS 

DNl=QV*SA*HA/( AMASS) 

DN2=G-HA/ ( VAIR**2 ) 

DN3=0. 

DN4=C0RD/ C2.*HA) 

DN5 = 0V*SA*C0RD*HA**2/ ( Y YI*32 . 2) 

DN6P=HA*G/(AMASS*VAIR**2) 

DN7P=ALT*(HA**2 )/((VAIR**2)*YYl) 

CALL INIT(T) 

CALL EOUIL 
PRMTC1 )=T 

CALL RKGS(PRMT, NDIM, IHLF, FCT ,OUTP,AUX) 

222 CONTINUE 
RETURN 
END 


Common Blocks 

This program has five common blocks. COMMON/UNKNOWN/ contains 
the solutions of the state variables at time step T. COMMON/VARIBL/ 
contains derivatives of the state variables at time step T. COMMON/ 
WINDS/ contains wind data at time step T. COMMON/TT/ contains wind 
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data. Unlabeled commons contain the aircraft data. The common blocks 
are used to transfer data from one subroutine to another. 


Common Block from Main Progra m 

COMMON CL,CLO,CLA,CLD£L,CLQ,CLADOT,CLG£ 

COMMON CD,CDO,CDA,CDA2,CDGE 

COMMON CM,CMO,CMA, CMDEL, CMQ , CM ADOT , CMGE 

COMMON QV, AMASS,G,PI , S A, CORD, ALT , HA, YYI, VAIR 

COMMON DN1,DN2,DN3,DN4,DN5,DN6,DN7,DN6P,DN7P 

COMMON FT,DELT,DEL, DELE, VA, GAMP 

COMMON TEST, NS, NP, NCO 

COMMON/UNKOWN/Y 

CDMMON/VARIBL/DERY 

COMMON/INI/VO 

COMMON/ WI NDS/ WX, WZ, WXT, rfZT,WXX,WXZ, WZX, WZZ, WXDOT, WZDOT , ZO , USTAR 
C0MMON/TT/AXAC41 , 11 , 2) , 0X,DZ 


Common Block from Subroutine IN IT 


COMMON CL,CLO,CLA,CLDEL,CLQ,CLADQT,CLGE 
COMMON CU,CDO,CDA,CDA2,CDGE 
COMMON CM,CMO,CMA,CMDEL,CMQ,CMADOT,CMGE 
COMMON Q V, AM ASS, G , PI, S A, CORD, ALT, HA, YYI, VAIR 
COMMON DN1,DN2,DN3,DN4, DN 5 , DN 6 , DN 7 , DN 6P , DN7 P 
COMMON FT, DELT, DEL, DELE, VA, GAMP 
COMMON TEST, NS, NP, NCO 
COMMON /UNKOWN/V ,GAM , 0 , AP, X,Z 

COMMON /VARIBL/VDOT, GMDOT,QDOT, APDDT, XDOT,ZDOI 
COMMON/INI/VO 

COMMON/WINDS/WX,WZ, WXT, WZT,WXX,WXZ,WZX, WZZ,WXDOT, WZDOT, ZO, USTAR 
COMKON/TT/AXA(41 , 11 ,2) ,DX,DZ 


Subroutine INIT 

Subroutine INIT is called by the main program. This subroutine 
defines and also initializes the various aircraft parameters. 

Nomenclature 

Subroutine INIT: 

APDOT a' Time derivative of angle of attack 
GAMD y Flight path angle in degrees 
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GAM 

Y 

Flight path angle in radians 

T 

t 

Time 

Q 

q 

Pitching rate 

FT 

f t 

Thrust 

HAD 


Reference altitude in meters 

Z 

z 

Vertical component of aircraft position, positive downward 

X 

X 

Horizontal component of aircraft position 

VA 

V a 

Aircraft velocity relative to air 

V 

V 

Aircraft velocity relative to ground 

VADOT 

• 

V a 

Time derivative of air velocity 

DEL 

6 

Angle between V a and V 

a 

GAMP 

Y 1 

Flight path angle in wind coordinate system 

VO 


Initial velocity 

Listing 

of 

Subroutine I N IT 


SUBROUTINE INITCT) 

COMMON CL,CLO,CLA, CL.DEL , CLQ , CLADOT , CLUE 
COMMON CD,CDO,CDA,CDA2,CDGE 
COMMON CM,CMO,CMA, C M D EL , C MQ , CMADOT , CMGE 
COMMON QV,AMASS,G,PI ,SA,CORD ,AL1 , H A , YYI , VAIR 
COMMON DN 1 , ON 2 , D N 3 , DN 4 , DN5,DNb,DN7,DN6P,DN7P 
COMMON FT, DELT, DEL, DELE, VA, GAMP 
COMMON TEST, NS, NP, NCO 
COMMON /UNKOWN/V ,GAM, 0 , A P, X,Z 

COMMON /VARIBL/VDOT,GMDOr,QDOT, APDOT, XDOT,ZDOI 
COMMON /INI/VO 

COMMON/ W I NDS/WX ,WZ,WXr,WZT,WXX,WXZ,NZX,rtZZ,WXDOT,WZDOT,ZO,USTAR 
COMMON /IT/ AX A (41 , 11 ,2) ,DX,DZ 
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NWR = 6 
APDQT=0, 

GAMD=0 • 0 

GAM=GAMD*PI/180. 

1 = 0 . 0 
3 = 0 . 
f r=i . 

HAD=HA*0. 3048 

KCK= 1 

HAM=HAD 

Z = - 1000,0/HA 

X = 0. 0 

XGM = X 

ZM = -Z 

CALL, WIND (XGM, ZM, HAM, KCK, WX , WZ , WXX , wXZ r WZX , WZZ) 

CALL CONV(VAIR,HAD,X,Z,XGM,ZM) 

10 F0RMATC2X,8E12.5) 

WRITE (6, 10) WX,WZ,WXT,WZT,WXX,WXZ,WZX,WZZ 

VADOT=0.0 

VA=1 .0 

V = WX* COS C GAM ) -WZ* SIN(GA.M) + ( (WZ* SIN(GAM)-wX* C0S(GAM))**2 
l+^A**2-CwX**2+WZ**2))**0.5 
DEL= ARSIN((WX* SIN(GAM)+WZ* COS CGAM) )/VA) 

GAMP=GAM+UEL 

VO=V 

WRITECNWR , 1 ) GAMD 

1 FDRMATC2X, • INITIAL FLIGHT PATH A NGLE= * , El 2 . 5 , ' DEG . ' , / ) 
RETUKN 
END 


Subroutine EQUIL 

Called from main program, this subroutine calculates initial 
values of the angle of attack, elevator angle, and thrust from the 
equations of motion by assuming equilibrium conditions: 


A(i,j) = B(i) , 


A 11 a ' + A 12 6 e + A 13 F T ~ B 1 * 


A 2i ot "F A 22 d £ "F ^2.3 T B 2 * 


A 31 0(1 + A 32 6 E + A 33 F T ~ B 3 5 
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from the following equations 


" D 1 V a C D “ D 2 sin + D 6 cos ( 6 T + a')F T - W = 0 , 

2 — 

D-j ^a “ ^2 cos + ^6 sin ( tS T + a '^T + W = 0 , 

°5 V a C m + °7 F I = 0 ’ 

where 

• • 

W = - W cos y' + W sin y' , 

A Z 

— • • 

W = W sin y' + W cos y' » 

A Z 

and the initial conditions, 

* • • 

Y = V = q = ot' ~ q = 0 . 

a 

Listing of Subroutine EQUIL 


SUBROUTINE EGUIL 

COMMON CL, CLO,CLA,CLDEL, CLQ , CLADOT , CLGE 

COMMON CD,COO,CDA,CDA 2 , CDSE 

COMMON CM, CMO, CM A, CM DEL, CMQ , CM ADOT , CMGE 

COMMON OV, AMASS,G,PI , SA , CORD , ALT , H A , 5 f Jf 1 , V AIR 

COMMON DNl,DN 2 ,DN 3 ,DN 4 ,DN 5 ,DNb,DN 7 ,DN 6 P,DN 7 P 

COMMON FT,DELT, DEL, DELE, VA, GAMP 

COMMON TEST , NS , N P , NCO 

COMMDN/UNKOWN/V , G A M , 0 ,AP,X,Z 

common/winds/wx , wz,wxr,<vzr,wxx,wxz,wzx,wzz, wXDOT, wzdot, zo, ustar 

DIMENSION A( 3 , 3 ),B( 3 ) 

NWR = 6 
AP= 0 . 

1 AP 1 =AP 
VA 2 =VA **2 

A(l, 1 )=-DN 1 *VA 2 *CCDA + CDA 2 *AP) 

A ( 1 , 2) =0 . 

A( 1 , 3 D=DN 6 P* COSCDELT+AP) 

A( 2 ,l)= DN 1 *VA 2 *CLA 
A ( 2 , 2 ) = DN 1 *VA 2 *CLDEL 
A C 2 , 3 ) -Dft 6 P* SIN ( DELT+A?) 

AC 3 , I )=Dub*VA 2 *CMA 
A( 3 # 2 )=DN 5 *VA 2 *CMDEL 
A C 3 , 3 ) =DN 7 P 
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WZDOT=WZT-V *(WZX*CDS(G4M ) -WZZ*SI N (GAM )) 

WXDOT^wXT+V *(WXX*CDS(G4M ) - WXZ*S1 N ( G AM )) 

Wl=WXDOT* CQS(GAMP) -WZDDT* SI N ( G A MP D 
B(1)=W1+DN1*VA2*CD0+DN2* SIN (GAMP) 

W2=rfXDOT* SIN(GAMP)+WZDDT* COS(GAMP) 

B(2)=-W2-DNl*VA2*CLO+DN2* COS (GAMP) 

B(3)=-DN5*VA2*CMO 
D1=B(3)/A(3,1)-B(2)/A(2, 1) 

C1=AC2,2)/A(2,1)“A(3,2)/A(3,1) 

C2=A(2,3)/A(2,1)-A(3,3)/A(3, 1) 

D2=B(3)/A(3, l)+A(3 f 2)*Dl/{AC3,l)»Cn 
C3 = A(3#2)*C2/(A(3, 1)*C1)-A13 , 3 ) / A ( 3 , 1) 

FTD=( A(1 , 1 )*D2-A(1 , 2)*01/Cl-BCn )/( AC1 ,2)*C2/C1-A(1 , 1 )*C3-A(1 , 3) ) 
AP=D2+C3*FTD 
DELE=-Dl/Cl-C2/Cl*FrD 
IF ( ABS(APl-AP) .GT.O.OOOOOnGO TO 1 
81 AP1=180. ♦AP/PI 

WRITE(NWK,2)AP1 , DELE , FTD 

2 F3RMATC2X, 1 EQUIL. ANGLE OF ATTACK= ' , El 2 . 5 , » EQUI L . ELEVAIOR AN GLE= 
1 , El 2 . 5 , 'EQUIL. THRUSr= ' ,E12.S,/) 

D.N6=DN6P*FTD 

DN7=DN7P*fTD 

WRITECNRR, 3 ) DN1,DN2,DN3 / DN4,DN5,DN6,DN7 

3 FORMAT (6X, 'DN1=',F12.5,/,6X,'DN2=',F12.5,/,6X,»DN3=',F12.5,/,5X,' 
1 D M4= > ,F12.5,/,6X, »DN5=' , F12.b,/,bX, 'DN6=* ,F12.5,/,bX, »DN7=' ,F12.5 
2 /) 

RETURN 

END 


Subroutine RKGS 

Subroutine RKGS evaluates a system of first order ordinary 
differential equations with given initial values by means of fourth 
order Runge-Kutta formulae in the modification due to Gill [24]. 
Accuracy is tested comparing the results of the procedure with single 
and double increment. This subroutine automatically adjusts the 
increment during the whole computation by halving and doubling. If 
more than ten bisections of the increment are necessary to get satis- 
factory accuracy, the subroutine returns with error message IHLF=11 
into main program [21]. (See pages 85 and 86 of this study for 

description of input/output parameters for this subroutine.) 
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Subroutine FCT 


Subroutine FCT is called by the subroutine RKGS. Subroutine 
FCT calculates the derivatives of the state equations at each time 
step. 

Listing of Subroutine FCT 


SUBROUTINE FCT ( T ) 

DIMENSION Y(6) , DER Y ( 6 ) 

COMMON CL,CLO,CLA,CLDEL,CLQ,CLADOT,CLGE 
COMMON CD,CDO,CDA,CDA2, CD GE 
COMMON CM,CMO,CMA,CMDEL,CMG,CMADOT,CMG£ 

COMMON QV,AMASS,G,PI,SA,C0RD,ALT,HA,YY1,VA1R 
COMMON DN1 , DN2 , DN3 , DN4 , DN5 , DNb ,DN7 , DN6P, DN7P 
COMMON FX,D£LT, DEL, DELE, VA, GAMP 
COMMON TEST, NS, NP, NCO 
COMMON/UN KUWN/Y 
COMMON/ VAK I BL/DERY 

COMMON/ WINDS/WX , WZ , WX T , *IZ I , WXX , taXZ , WZX , WZZ , rtX DOT , WZDUT , ZO , U ST AR 
E'JU I VALENCE CDERYC 1 ) , V DO T ) , C DER Y C 2 ) , GMDOT ) , CDERYC3) , ODOT) , 
1(DERY(4) , APDOT) , ( DER Y C 5 ) , XDOT) , (DERYtb) , ZDOT) 

EOUIVALENCE(YCl) , V) , C Y C 2 0 , GAM) , (YC3),Q),(YC4),AP),(Y(5),X), 
1(Y(6),Z) 

CALL ARCOEF 

DERY C 1 ) = VDOT , DERY t 2 ) = GMDOT , DER Y ( 3 ) =0 DOT, DERY C 4 ) =APDOT , DEKY ( 5) 
XDOT, DERY (6 )=ZDOT 

Y(1)=V ,Y(2)=GAM ,Y(3)=3, Y ( 4 ) =AP , Y ( 5 ) =X , Y ( 6 ) =Z 
101 DEL= ARS1N((WX* SI N ( G A M ) + W2* COS ( G A M ) ) / V A ) 

G AMP= Y ( 2 ) + DEL 
DERY ( 5) =V* COS C GAM ) 

DERY ( 6 ) s-V ♦ SIN ( GAM ) 

VA=( (DERY(5)-WX)**2+(DERY(6)-wZ)**2)**0.5 
WZDOTs WZT- V *(WZX*COS(GAM ) -lfiZZ*SI N ( GAM )) 

WXDDT=HXT+V * (WXX*COS(GAM ) - ft XZ*S1 N ( GAM )) 

DERY (1)=-DN1*(CD*CCJ$(DEL)+CL*S1N(DEL))*VA**2-DN2*SIN(Y(2))+DN6 
l*COS(DELT+Y (4)+DEL) 

DERYC2)=DN1*VA**2*(CL*CDS(DEL)-CD*SIN(DEL) ) /Y ( 1 ) -DN2*COS ( Y C 2 ) ) / 
1Y(1 )+DN6*SIN(DELT + Y(4) + 0EL)/.YCl ) 

DERY (3)=DN7*FT+DN5*VA**2*CM 

GMPDOT =DN1*CL*(VA)-DN2* COS(GAMP)/VA +DN6*FT* 5IN (DELT+Y C 4 ) ) / 
1 V A +CWXDOT* SIN(GAMP)+WZDOT* COS t GAMP )) /VA 
DERY(4)=YC3)- GMPDOT 
10 RETURN 
END 
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Subroutine ARCOEF 


Subroutine ARCOEF is called by subroutine FCT. Subroutine 
:al 

time step: 


ARCOEF calculates the aerodynamic coefficients CL, CD, and C at each 

m 


C^ — C L + C L ex 1 + C, 6,- + ot?" C 


a 


+ r + a r 

L x ~E ' 2V 2V A L« r ’ 


a q a a 


'GE 


C n - C n + C n a 1 + C n 0 a 1 + AC n , 


'D D D 
o a 


D 2 

a 


GE 


C=C + C a 1 + C 6,- + c + 
m m m ~ 


Ca 


a 


\J | — I At • Vj 1 Oil C /\ C 

m, E 2V m 2 V_ m- m 


[19] 


a q a a 


'GE 


For the augmentor-wing STOL aircraft, the equations of the 
aerodynamic coefficients vary significantly from those presented above, 
and are expressed as: 


+ £S- C, . , 


C L + + a + 2V 2 V L 

L l WB l H l H L q ^ v a a 

°° o a 


C D C D + C D + &C D r| - ’ 
o a GE 


Z W 

C = C + C a 1 + C. ( — cos a' + — sin a' ) 
m m o m a L WB c c 


£ Z 

C n ( sin a 1 - - 3 - cos a 1 ) + ( C L + C. a' + C L 6 ^ ) 
c c H o H a 6 E 


sin a' - ^ cos a' ) + C m + ff- C m . + AC m _ , 
c c a q a a 


'GE 
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where 


C, = C, + C. a' , 

l wb L m n l wb 

00 o a 

prior to inclusion of ground effect terms [19]. 

Since the ground effect, terms are negligible for conventional 
aircraft (see Figure 23, page 44), they are set equal to zero for 
convenience. Their inclusion is a trivial task when given the par- 
ticular equations for each of the terms, as presented in Chapter IV. 

Listing of Subroutine ARCOEF 


subroutine: arcoef 

COMMON CL,CLO,CLA,CLDELi,CLiQ,CLADUT,CLGE 
COMMON CD,CDO,CDA,CDA2,COGE 
COMMON CM,CMO,CMA,CV!DEL,CMQ,CMADOT, CMGE 
COMMON GV, AMASS, G, PI ,SA, CORD ,ALT, HA , Y YI , VAIR 
COMMON DN1 ,DN2,DN3, DN4, DN5,DN6, DN7, DN6P,DN7P 
COMMON FT, BELT, DEL/DEL£,VA, GAMP 
COMMON TEST,NS,NP,NCQ 
COMMON /UN KOWN/V, GAM,Q, A?,X,Z 

COMMON /VARI6L*/VDOr,G v iD3T,QDOT,APDOr,XL)QT,ZDOT 

CLGE=0. 

CDGE=0. 

CMGE=0. 

CU=CLO+CLA*AP+CLDEL*DELE+DN4*(CLQ*G+CLADOr*APDOT)/VA+CLGE 

CD=CDO+CDA*AP+CDA2*AP**2+CDGE 

CM=CMO+CMA*AP+CMDEL*DELE+DN4* ( CMQ*Q+CMADQT* APDOT) /VA + CMGb 

RETURN 

EftiD 


Sub routine OUTP 

Subroutine OUTP is called by the subroutine RKGS at every time 

step. Subroutine OUTP is an external output subroutine used to print 

the results. Subroutine OUTP in turn calls subroutine WIND and updates 

the WIND speed parameters. Also, if the automatic control system is 

desired, subroutine OUTP calls subroutines MODE and GAINS which 
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computes the control parameters Fy and 6 £ . Subroutine OUTP terminates 
the integrating subroutine RKGS when the aircraft reaches ground, 
z = 0, by changing the value of PRMT(5) to a non-zero value. 

Listing of Subroutine OUTP 


SUBROUTINE OUTP ( T , I HLF , N DI M , PRMT ) 

DIMENSION Y ( 6 T -PRMTC5) . DER Y ( 6 ) 

COMMON CL,CLO,CLA,CLDEL,CLQ,CLADOTrCLGE 

COMMON cd,cdu,cda,cda2 , cdge 

COMMON CM ,CMO , CMA , CM DEL, CMQ , CMADOT , CMGE 

COMMON QV , AMASS, G, PI , S A , CORD , ALT , H A , Y t I , V AI R 

COMMUN ONI ,DN2,DN3,DN4,DN5,DN6,DN7 , DN6P,DN7P 

COMMON FT,DELT, DEL, DELE, VA, GAMP 

COMMON T£ST,NS,N2,NCO 

COMMON/UNKOWN/Y 

COMMON /VARIBL/DERY 

COMMON/INX/VO 

COMMON/WINDS/WX , WZ , WXT , WZT , wXX , WX Z , NZX , WZZ,WXDOT,WZDOT,ZO,USTAR 
COMMON/ TT/AXAC 41 , 1 1 ,2) , DX ,DZ 

EOUI VALENCE ( DERY (1 ),VD0T),(DERY(2) , GMDOT ) , ( DER Y ( 3 ) ,QDOT) , 
1CDERYC4) , APDOT) , ( DERY ( 5 ) , XDOT ) , ( DERY ( 6 ) , ZDOT) 

EQUIVALENCES 1) ,V),CYC2),GAM),CYC3),Q),CYC4),AP),CYCS),X), 

1 (Y(6) ,Z) 

C 

33 NNR = 6 

HAD=HA*0. 3048 

XGM=X*HAD 

ZM=-Z*HAD 

1FCN2.GT.1) GO TO 4 

IIC = 0 

M1=N2 

M2 = N2 

M3 = N2 

M4 = N2 

MS=N2 

M6 = N2 

ICHK=0 

N 1 = 2 0 

M7 = 0 

r r= o . 

4 N 2 = 2 
DT=T-TR 
IR = T 

IFCNCO.LT. 10) GO TO 100 
NCO = 0 

WRiTECNWR, 12)T 

12 FORMATC5X, 'TIME = *,F12.6) 

WRIIE{NWR,10)(Y(I),I=1,NDIM) 

10 FORMAT ( 2X , ' V = 1 , F9 . 5 , 2X , ' GAMMA = ' , F9 . 5 , 2X f ' PITCH RATE 

1 ' f F9 . b ,/,2X, 'ALPHA PRIME = ' , t 9 . 5 , 2X , 'HORIZONTAL POSITION = 
2F9.b,2X, ' VERTICAL POSITION = '.F9.51 
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WRITE(NWK, 13) (DERV (I) , 1 = 1 , NDIM) 

13 FORMAT ( 2X , ' VDOT = ',E15.5,2X,' GhDDT = 1 ,E15.5,2X, * QDOT = ',E15 .5 
1,/,2X, ' AFOOT = ',E15.5 ,2X,»XDGT = * , El 5 . 5 , 2X , ' ZDOT = ',E15.b) 
WRITE (NWK, 30) WX,WZ>WXDDT., W2D0T,GAMP, VA, DEL 
30 F0RMAT(2X, ' WX = ' , El 2 . 5 , 1 X , ' W2 = • , El 2 . b , 1 X , » WXDDT = • , El 2 . 5 , 1 X , 

2'WZDOT =',E12.5,1X,'GAM?=',E11.4,1X, ' VA= » , Ell. 4, IX, 'DEL =' ,E11.4) 
WSITEC6 , 81 ) CL, CD, CM 

81 FORMAT (6X, ' CL= ' ,E12.5, ' CD=' ,E12.5, 'CM= • ,E12.5) 

WRITE(6, 1001) XGM , ZM 

1001 FORMAT C2X, «XGM= « ,E12. 5, ' ZM = ' ,E12.5) 

WRITE (6,8888) FT , DELE , THET 

8888 FORMAT ( 4 X , ' FT= ' , El 2 . 5 , 2 X , < DELE= ' , El 2 . 5 , 2X , ' THET= • , El 2 . 5 ) 

100 NC0=NC0+1 
b CONTINUE 
X3M = X 
ZM = -Z 
KCK = 2 
HAM=HAO 

CALL MODE ( M 1 , M 2 , M 3 , M4 , Z , ZDOT , X , VO , V , G AM , DT , THET ) 

CALL GAINS (THET, Mb, M6) 

CALL WIND (XGM, ZM , HAM, KCK, wX , WZ, WXX, WXZ, rtZX, WZZ) 

CALL CONV (VAIR, HAD, X,Z,XGM,ZM) 

CHK1=X*HAD 
CHK2=-Z*HAD 
ENDZ=10.*DZ 
ENDX=40. *DX 

IF(CHK1 .GI.ENDX) GO TO 300 
1F(CHK2.GT.ENDZ) GO TO 300 
1 F ( Z ) 200,300,300 
300 PRMT ( 5 ) = 1 • 0 
200 CONTINUE 

IF(IHLF-IO) 1 ,1,2 
2 WRITECNWK, 1 1 )IHLF 

11 FQRMAT(2X, 'ERROR IHLF = ',14) 

1 CONTINUE 

RETURN 
END 


Subroutine MOD E 

Subroutine MODE is called by the subroutine OUTP. The input 
parameters , M2, M^, and M^ are initially zero. Other input param- 
eters are Z, ZDOT, X, VO, V, GAM, and DT. The output parameter is 
THET. Subroutine MODE in turn calls one of the subroutines HOLD, 
CAPTR, TRACK, or FLARE, depending upon the position of the aircraft 
mode control logic. 
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Listing of Subroutine MODE 


SUBROUTINE MODE ('ll ,M2 f M3,M4,Z f ZD0T,X r VD,V, GAM, DT, THET) 

IF(ABSCZ) . LE.0.2) GO TO 23 

IFCX.LT.3.0) GO TO 20 

IF(ABSCGAM) .GE. 0.047 ) GO TO 22 

IFCX.GE.3.4) GO TO 22 

IFCX.GE.3.0) GO TO 21 

20 CALL HOLDCMl ,Z, DT,TH£T) 

GO TO 24 

21 CALL CAPTR(M2,Z,ZDOr,DT,X,GAM,VO,THET) 

22 CALL TRACKCM3 , Z, ZDDT, DT , X , GAM , VO , THET ) 

GO TO 24 

23 CALL FLARE(M4,Z,ZD0T,V,DT,THET,1C) 

IFCIC.EQ.I) GO TO 22 

24 CONTINUE 
RETURN 
END 


Subroutine HOLD 

Subroutine HOLD is called by the subroutine MODE. The input 
parameters are I, Z, and DT. Input parameter I is initially zero and 
then set to I = 2 in this subroutine. The output parameter is THET. 
The following difference equations are solved in this sub- 

routi ne: 


AZ = 1 - HR s 

a. DT 

C 1 = e 1 , 

C 2 


C 3 « K 2 (a 2 DT - 1) , 


™j " C 1 TH j-l + C 2 AZ j-l ’ 


THET. = THET j _ 1 + K 2 THj + C 3 TH j _ ] . 
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Listing of Subroutine HOLD 


SUBROUTINE HOLDCI,Z,DT,THET) 
C HOLD MODE 

HR = Z 
T1 A = 1 • 

T1K=0 . 1 
T2A=0 , 1 
T2K=0.05 

IF ( 1 , GT . 1 ) GO TO 5 
T1R=0. 

T2R=0. 

T3H=0. 

1 = 1 + 2 

5 CONTINUE 
T1=Z-HR 

T1C=EXP(-T1A*DT) 
riD=TlK*(l.-TlC)/Tl A 
T2=(T1C*T2R+T1D*T1R) 

T2C=T2K* ( T2 A*DT- 1 . ) 

T3=T3R+T2K*T2+T2C*T2R 

THET=T3 

T3R=T3 

T2R = T2 

T1R=T1 

return 

END 


Subroutine CAPTR 

Subroutine CAPTR is called by subroutine MODE. Input param- 
eters for this subroutine are I, Z, ZDOT, DT, and VO. Input parameter 
I is initially zero and then set to I = 2. In this subroutine the 
output parameter is THET. 

The following difference equations are solved in this sub- 
routine: 

™j = C 1 TH j-l " C 2 TH j-2 + C 3 AZD j-l + C 4 ZD j-2 + ATHETP ■ 

AZD. = HCD - ZDOT , 

J 

and where 
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c 

c, 

c 

c 


= 1 + e^ , 

- e- DT / T , 

= K 3 k 4 [dt + T (C 2 - 1) ] , 
= K 3 K 4 [DT - C 2 (T + 1) ] . 


Listing of Subroutine CAPTR 


SUBROUTINE CAPTR(I,Z,ZDDT,DT,X,GAM, VO,THET) 

T AO= 0 ,01 
CT1K=0. 1 
THETP=-0 . 05 
I E ( I . GT , 1 ) GO TO 5 
CT1R=0. 

CTlRRrO. 

CT2R=0. 

CT2RK=0. 

1 = 1 + 2 

5 CONTINUE 

HCD=VQ*SIN(4. 71239E-2) 

CI1=HCD-ZD0T 
CT1C1=1 .+EXPC-DT/TAO) 

CriC2=EXP(-DT/TA0J 
CIlDlrCTlK*CDT + TA0*CCTlC2-l . ) ) 
CTlD2=CTlK*(TAO-CTlC2*(Dr+TAO) ) 

Cr2=CTlCl*CT2R-CTlC2*Cr2RR+CTlDl*CIlR+CTlD2*CTlRR 

CT2=CT2*U.0l 

THET=CT2 

CllRsCTl 

CT1RR=CT1R 

CT2R=CT2 

CT2KR=CT2R 

THET=THET+THETP 

RETURN 

END 


Subroutine TRACK 

Subroutine TRACK is called by subroutine MODE. Input param- 
eters for this subroutine are I, Z, R, DOT, X, GAM, and VO. Input 
parameter I is initially zero and then set to I = 2 in this subroutine. 
The output parameter is THET. The following difference equations are 


solved in the subroutine: 
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e. = GAM. - GAMR , 

J J 


TH 


j = K 4 (Zj - ZR) + Kg(ZDOT. - ZD) , 


TT1 j C 1 TT1 j-l + C 2 £ j-l ’ 


TT2 j “ TT2 j-l + K 6 TT1 + C 3 TT1 j-! * 


THET. = TH. + TT2 . , 

J J j 


where 

c i 


= e 




= rO-Cj , 


C 3 = M a 3 DT ~ 


Listing of Subroutine TRACK 


SUBROUTINE TRACKCI , Z, ZDOT ,DT,X, GAM, VO, THET) 
C TRACKING MODE 

TT3K=0.01 
TT3K=0 .005 
TT3K=0.0b 
TT4K=0 .05 
TT4K = 1 . 

TT4K=0. 00325 
I AU = 0 , 1 
TAU=0.01 

rriA=i ,/tau 
,/tau 

TT2K=1 . 

TT2K=0. 125 

IT2K=0. 00125 

TT2 A=0 .01 

GAMR=-0. 0471239 

HR=-(63-X)*TAN(-GAMR) 

DFH=Z-HR 
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IF ( I . GT . 1 ) GO TO 5 

rr2R=0. 

TT1R=0. 

EPSK=0. 

1 = 1 + 2 

5 CONTINUE 

DRH=ZDOT-VO*SIN(-GAMR) 

EPS=G AMR-GAM 
rTlC=EXPC-TTl A+DT) 

TT1D = T'11K*( l.-TTlC)/rriA 

TI1=TT1C*TT1R+TT1D*EPSR 

TT2C=TT2K»(TT2A*DT-1 . ) 

rr2 = TT2R + TT2K*TTl + Tr2C*m R 

THETT=TT2+TT3K*DFH+rT4K*DSH 

THET=TH£TT 

TT1R=TT1 

EPSR=EPS 

I T2R=TT2 

C RATE LIMITER 

IF(GAMR-GAM)! ,4,2 

1 IF (THET ) 4 , 4 , 3 

2 IF ( THET ) 3 , 4 , 4 

3 XHET=-0. S+THET 

4 RETURN 
END 


Subroutine FLARE 

Subroutine FLARE is called by subroutine MODE. Input param 
eters for this subroutine are M, Z, ZDOT, V, DT, and IC. Input 
parameters IC and M are initially zero and then set to IC = 1 and 
M = 2 in this subroutine. The output parameter is THET. In this 
subroutine initially when M = 0, three different reference heights 
(HSTP, HF, and HRMP) for FLARE are calculated. Depending upon the 
altitude of the aircraft Z, the output parameter THET is calculated 
The following difference equations are solved in this subroutine: 


P5 j 

= (PSj., 

PI . 
J 

' (ZD0T j 

P2 j 

■ P1 .T K 1 


- THRj_ 1 ) DT , 

- ZD ) AA2 + Z . 

vJ 

AA , 


9 
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P3 j ■ P3 j-1 + P2 j ’ 

P 6j C] + KfTHPj - THPj_-|) , 

THET j ' P2 j - P3 j + P5 j + P6 j • 
where 

r -A6*DT 

Ci = e 


List ing of Subrout ine FLARE 


SUBROUTINE FLARECMfZ/ZDOT/V ,DT,TH£T,IC) 

C FUAHE MODE 

IFCM.GT. 1 ) GO TO 2 
TCL. = 0. 

P4=0.0 
P5=0.0 
Pfa = 0.0 
P3R=0. 

P5R=0. 

P6K=0. 

TriRR=0. 

thpk=o. 

HSTP=3 . 5*ZDOT/V 
HRMP=0.6*HSTP 
HF=0 . 8*HSTP 
THP=0.9*ZDOT/V 

THK=C(0.255*Z DOT/V)-0. 0001315)* DT 
WKITE(6,12)Z,HSTP,THP,THR 
12 F3RMATC/ , 5X, 4E20.6, /) 

2 M=2 

zoorc=o.oi 
AK1= 0.00024 
AA1=0. 00333 
AA2=5.0 
AKP=2. 

A AP=5 . 

IF(ABSCZ) .GT.HSTP) GO TO 6 
IF(ABSCZ) .GT.HF) GO TO 5 
IF( ABS(Z) .GT.HRMP) GO TO 4 
P5=(P5R+TriRR)*0.01 
P5R=P5 
TrlRR = THR 
TCL=TCL+DT 
4 CONTINUE 
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P1=(ZD0T-ZD0TC)*AA2+Z 
I F ( TCL • GT .4.0) PS = 0. 

PI 1 = P 1 * AK 1 
IFCPU.GE.0.)AA = AA1 
IF'CPll . LT . 0 . )AA = AAl/2, 

P2=P11*AA 
XMR=P2 
P3=P3R+XNR 
P3R = P 3 
P4=P2-P3 

b CONTINUE 

ClP=EXPt-AAP*DT) 

P6=ClP*PbR-*-AKP*(THP-THPR) 

PbR=P6 
THPR=THP 
THET=P4+P5+P6 
C RATE LIMITER 

IFCTHET.GT.O. 1 ) T HE T = 0 . 1 

1F-C THET.Lr. 0.2E-04I rHET = 0. 2E-04 

GO TO 7 

6 IC = 1 

GO TO 8 

7 10 = 0 

8 CONTINUE 
RETURN 
END 


Subroutine GAINS 


Subroutine GAINS is called by subroutine OUTP. Input param- 
eters for this subroutine are THET, Ml, and M2. Aerodynamic 
coefficients, AP, APDOT, DEL, Q, and V A are brought in by common 
statements. This subroutine in turn calls subroutines SERI and SER2 
calculates elevator and thrust signals. The following equations are 
solved in this subroutine: 


F tc = KT1 VA - KT2 Z/V + KT3 X/V + K4 THET , 

6 ec = KD1 VA - KD2 Z/V - KD3 X/V - KD4 THET , 

where 
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and 


KT1 


G io G n ~ G i 2 


fH ^ 
n 9 

KT2 


D2 0 0 


H 1n 


z: 



10 

KT3 


0 -D2 0 


H n 

KT4 


0-10 


. 


kdi' 


- G 10 G 11 G 12 


V 

KD2 


-D2 0 0 


H 7 

KD3 


0 -D2 0 


H 8 

KD4 


. 0 -1 0 . 


l 


Listing of Subroutine GAINS 


SUBROUTINE GAINS(IHET,M1 , M2 J 

COMMON CL,CL.O,CLA,CLDEL,CLQ, CLADOT,CL.GE 

COMMON CD,CDO,CDA,CDA2, CDGE 

COMMON CM,CMO,CMA,CMDEL,CMQ,CMADOT,CMGE 

COMMON Q V, AM ASS, G, P I, S A, CORD, ALT, HA, Y ^ , VA1R 

COMMON DN1,DN2,DN3,DN4,DN5,DN6,UN7,DN6P,DN7P 

COMMON FT,DELT,DELi, DELE, V A , GAMP 

COMMON TEST, NS, NP, NCO 

COMMON/UnKOWN/V ,GAM ,3,AP,X,Z 

COMMON /VARIBL/V DOT , GM DOT , QDOT , APDOT , XDOT , ZDOT 

COMMON/GNS/AFK1 , AFK2 , AFK3 , AFK 4 , ADK1 , ADK2 , ADK3 , ADK4 

F 1 =AP 

F2=APOOI 

F3 = D£Li 

F 4 = V A 

F5 = 0 

FA=COSCF3) 

FB=SINCF3) 

FC=COS(DELT+AP+DEL) 

FD=SINCDELT+AP+DEL) 

C1=DN1*CDA 

C2 = DN 1 *CDA2 

C3=DN1 *CLA 

C4 = DN 1 *CLDEL 

C5=ON5*CMA 

C6=DN5*CMDEL 

C7=DNl*CDO 

C8=DN1 *CLO 

C9 = DN1*DN4*CL.U 

C10=DN1*DN4*CLAD0T 

Cl l=DNS*CMO 

C12=DN5*DN4*CMQ 

C13=DNS*DN4*CMADOT 
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G1=DN7 

G2=~C1*FA“C2*FA*F1"C3*F3 

G3=-C4*FB 

G4=DN6*FC 

GS=C3*FA-C1*FB-C2*FB*F1 

G6=C4*F A 

G7=DN6*FD 

G8 = Cb 

G9 = C6 

G10=C7*F4*FA+C8*F4*FB+C9*F5*FB+C10#F2*FB 

G1 1=C8*F4*FA+C9*F5*FA+C1 0*F2*FA-C7*F4*FB 

G12=C1 1*F4+C12*F5+C1 3* F2 

H1=G6*G1-G7*G9 

H2=G5*G1-G7*G8 

H3=G5*G9-G6*G8 

H4=G2*H1-G3*H2+G4*R3 

Hb=F4*H4 

H6=H2/(F4*H5) 

H7=(G2*G1-G4*G8) /(F4*H5) 
H8=(G2*G7-G4*G5)/(F4*H5) 

H9=H3/H4 

H1Q=(G2*G9-G3*G8)/H4 

HU = (G2*Gb-G3*G5)/H4 

ADK1=G12*H8-G11*H7-G10*H6 

ADK2=-DN2*Hb 

ADK3=-DN2*H7 

ADK4 = -H7 

AFK1=G10*H9+G11*H10-G12*H11 

AFK2=DN2*H9 

AFK3=-UN2*H10 

AFK4=-H10 

GA1=ADK1*VA-ADK2*ZD3T/V - A OK 3 *XDOT/ V - ADK4 * THET 
GA2=AFK1*VA-AFK2*ZDDT/V+AFK3*XDGT/V+AFK4*TH£T 

CALL SER1(M1,GA1,DELE,DT) 

CALL SER2(M2,GA2,Fr,Dr) 

RETURN 

END 


Subroutine SERI 


Subroutine SERI is called by subroutine GAINS. Input param- 
eters for this subroutine are M, XN, and DT. Input parameter M is 
initially zero and then set to M = 2 in this subroutine. Output 
parameter is control variable DELE. This subroutine calculates the 
control variable DELE from the control signal XN (calculated in sub- 
routine GAINS) using an elevator control servo. The difference 
equations solved in this subroutine are as follows: 
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T1 j = Ct T1 + K/t (XN. - XN...,) . 

DELE. = C 2 DELEj., - C 3 DELE.., + C 4 Tlj + C g T lj _ 1 + C fi Tl^ . 
where 

-a-,/T DT 
C ] = e 1 

-a^DT 

C 2 = 2 e cos b , 

-2a ? DT 

C 3 = e , 

C 4 = K 2 /RT (1/RT - sin ip/b) , 

C 5 = K 2 /RT (RD1 + sin(ijj DT) + RD2) , 

a 9 DT 

C 6 = K 2 /RT e L (RD3 - RD4) . 


Listing of Subroutine SERI 


SUBROUTINE SER 1 ( M , X N , DELE , DT ) 

C SERVO SYSTEM 

IF(M.GT.O) GO TO 2 
T 1 R = 0 • 

TlRrt=0. 

T2R=0. 

T2RR=0. 

XNR=0. 

2 M = 2 
AK1=1. 

AAI=0. 00036 
TA1=0.995 

or=o.oi 

AAT=AA1/TA1 

AKT1=AK1/TA1 

C1=EXP(-AAT*DT) 

AK2=390000. 

AA2=883. 17469 
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BB=1 .5707963 

RT=SQRT( (AA2*AA2)+(BB*BB) ) 

AKT2= AK2/RT 
PSI=ATAN(-SB/AA2) 

XXP=EXP(-AA2*DT) *COSC SB) 

Cr2=EXP(-2.*AA2*DT) 

Dn = AKT2*C(l./RT)-(SIN(PSI)/BB) ) 

CH = 2. *XXP 
RD1=-2.*XXP/RT 

R02 = EXP(-AA2*DT)*SIN CBB*Dr + PSI) 

Dr2 = AKX2*CRDl+5IN CPSIfOD +RD2) 

RD3=EXP(-AA2*DT)/RT 

RD4=-SINCBB*DT+PSI)/B8 

DI3=AKT2*EXP(-AA2*Dr)*tRD3-RD4) 

ri-Cl*TlH+AKTl*(XN-XNR) 

f2=CTl*T2K-CT2*T2RR+DTl*Tl+DT2*TlR+DT3*TlRR 

riRR=TlR 

T1 R=T 1 

T2RR=T2R 

T2R=T2 

X N R = X N 

DEBE=T2 

RETURN 

END 


Subroutine SER2 

Subroutine SER2 is called by subroutine GAINS. Input param- 
eters for this subroutine are M, XN, and DT. Input parameter M is 
initially zero and then set to M = 2 in this subroutine. Output 
parameter is control variable Fy. This subroutine calculates the 
control variable from the control signal XN (input from subroutine 
GAINS) using a thrust control servo. The difference equations solved 
in this subroutine are as follows: 


T1 j 

c i 

T1 M 

FT 0 

■ C 2 

FT J-1 

where 




+ K/t (XN. - XN j _ 1 ) > 

- C 3 FT j-2 + C 4 T1 j + C 5 T1 j-1 + C 6 T1 J-2 


= e 


-a-j/x DT 


= 2e 


-a2DT 


cos b 
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2a«DT 

C 3 e 

C 4 = K 2 /RT (1/RT - Sin ip/b) , 

C c = K 0 /RT (RD1 + sin DT) + RD2) , 
b c- 

-a 9 DT 

C 6 = K 2 /RT e (RD3 - RD4) . 


List ing of Subroutine SER2 


SUBROUTINE SER2(M,XN, FT,DT) 

C SERVO SYSTEM 

IF(M.GT.O) GO TO 2 
T1K=0. 

T1RK=0. 

T2R=0. 

X2RR=0 • 

XNR=0. 

2 M = 2 

DT = 0 . 0 1 
A K 1 = 1 . 

AA1=0. 00036 
TA1=0 . 995 
AAT=AA1/TA1 
AKT1 = AK1 /l'Al 
Cl=EXPt-AAT*DT) 

AK2=390000. 

AA2=883. 17469 
BB=1 .5707963 

RT = SORT( C AA2 + AA2 ) +(BB*BB) ) 

AKT2=AK2/KT 

PSI=ATANC-BB/AA2) 

XXP=EXPC-AA2*DT)*C0S(BB) 

CT2=EXP(-2.*AA2*DT) 

Dn = AKT2*C ( 1 ./RD-CSIMCPSI) /BB) ) 

CT1=2.*XXP 

RD1=-2.*XXP/RT 

RD2 = EXPOAA2*DT)*SlNCBB*Dr + PSI) 

0r2 = AKT2*(RDl + SIN(PSI^DD+RD2) 
RD3=EXP(~AA2*DT) /RT 
RD4=-SINCBB*DT+PSI J /BB 
DT3=AKT2*EXP(-AA2*Dr) *( RD3-RD4) 
T1=C1*T1R+AKT1*(XN-XNR) 

T2=CT1 *T2R-CT2*T2RR+DTl*Tl+DT2*TlR+Dr3*TlRR 

T1KR=T1R 

T1R=T1 

T2RR=T2R 

T2R=T2 

XmR=XN 

FT=T2 

RETURN 

END 
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2a,DT 

C 3 -e , 

C 4 = K 2 /RT (1/RT - sin ifi/b) , 

C 5 = K 2 /RT (RD1 + sin (\p DT) + RD2) , 
-a ? DT 

Cg = K 2 /RT e (RD3 - RD4) . 


Listing of Subroutine SER2 


SUBROUTINE SER2(M,X'J,FT,DT) 

C SERVO SYSTEM 

IF(M.GT.O) GO TO 2 
T1K=0. 

T1RR=0. 

T2R=0. 

T2RR=0. 

XNK=0, 

2 M = 2 

DT=0 . 0 1 
AK1=1 . 

A A 1 = 0 . OOQ 3 6 

TAl=0.9y5 

AAT=AA1/TA1 

AKT1=AK1/1'A1 

Cl=EXPl-AAT*DT) 

AK2=390000. 

AA2=b83. 1746y 
66=1.5707963 

RT=SQRT( (AA2*AA2)+(BB*BB) ) 

AKT2=AK2/KT 

PSI = ATAN C-BB/AA2) 

XXP=EXP(-AA2*DT)*C0S(BB) 

CT2=EXP(-2.*AA2*DT) 

DT1=AKT2*( (1./RT)-(S1MCPSI)/BB) ) 

CT1=2. *XXP 
RDl=-2 . +XXP/RT 

RD2=EXPC-AA2*UT) *S 1 N ( B B * D T+PS I ) 

0X2 = AKT2*(RD1+SIN(PSI*DTMRD2) 

R03=EXP(-AA2*DT)/RT 

RD4=-SINCBB*DT+PSI)/BB 

DT3 = AKT2*EXP(-AA2*DD *( RD3-RD4) 

n=Cl*TlR + AKTl*(XN-XNR) 

T2=CT1 *T2R-CT2*T2RR+DT1 *T1 +DT2 *T 1 R+D T 3*T 1 RR 

T1KR=T1R 

T1R=T1 

T2RR=T2R 

T2R=T2 

XNR=XN 

FT = T2 

RETURN 

END 
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